Semi-Dilute Dumbbells: Solutions of the Fokker–Planck Equation

Abstract

Spring bead models are commonly used in the constitutive equations for polymer melts. One such model based on kinetic theory—the finitely extensible nonlinear elastic dumbbell model incorporating a Peterlin closure approximation (FENE-P)—has previously been applied to study concentration-dependent anisotropy with the inclusion of a mean-field term to account for intermolecular forces in dilute polymer solutions for background profiles of weak shear and elongation. These investigations involved the solution of the Fokker–Planck equation incorporating a constitutive equation for the second moment. In this paper, we extend this analysis to include the effects of large background shear and elongation beyond the Hookean regime. Further, the constitutive equation is solved for the probability density function which permits the computation of any macroscopic variable, allowing direct comparison of the model predictions with molecular dynamics simulations. It was found that if the concentration effects at equilibrium are taken into account, the FENE-P model gives qualitatively the correct predictions, although the over-shoot in extension in comparison to the infinitely dilute case is significantly underpredicted.

1. Introduction

Spring bead models prove to be an effective way to construct constitutive equations for polymer melts. Perhaps the most widely used with success is the finitely extensible nonlinear elastic (FENE) dumbbell model, which exhibits a viscous response that arises from Brownian motion and an elastic response that arises from the spring connecting the beads [1]. However, single dumbbell models are just a simple approximation to polymer melt solutions and do not reflect all the complex physical processes that can occur. Augmented dumbbell models have been developed that include additional effects such as internal viscosities [2,3], additional hydrodynamic interactions [4,5,6], and aniostropic drag [7]. These advances resulted in a computationally viable constitutive model, allowing improved predictions and greater understanding of single polymer experiments [8,9,10,11].

The aim of this article is to further investigate the mean-field model proposed by Schneggenburger et al. [12], which also incorporated concentration effects and weak shear. This model is revisited in [13] to include steady uniaxial elongation within the Hookean regime. The behavior of the FENE dumbbell model without the mean-field term is well understood and has been extensively studied [14,15]. The mean-field model under uniform weak shear-flow conditions gives good predictions of the orientation angle and qualitative predictions of the radius of gyration when compared with light scattering experiments [12,16]. Concentration effects have been quantified experimentally under equilibrium conditions and for shear flows [16,17,18]. Concentration-induced effects were found to be more significant in elongational flows than in shear flows [19]. A study using a four-roll mill found that increases in concentration effects result in a reduction of the extension due to polymer entanglement [20]. A later study using a FENE-CR fluid to model the full four-roll mill system found qualitative agreement with experiments [21]. This reduction in extension due to semi-dilute behaviors has also been reported experimentally [22], in Brownian dynamical simulations [23] and from modified dumbbell models [24].

In this paper, we will extend the investigations in [12,13], the authors of which solved the constitutive equations for the second moment. We will include the effects of strong shear and large elongation. Further, we will explore the form of the probability density function, $\psi$, which permits the computation of any macroscopic variable, thus enabling direct comparisons of the FENE dumbbell model with molecular dynamics simulations. The paper is structured as follows. The governing equations are introduced in Section 2. The equations are solved Section 3, with results presented for the dumbbell model with mean-field force under strong shear and elongation. Finite element solutions of the FENE model incorporating the Peterlin closure approximation (FENE-P) are also presented, along with results of molecular dynamics simulations. Results are discussed and conclusions drawn in Section 4.

2. Materials and Methods

The Governing Equations

We will model the fluid as a suspension of dumbbells embedded in a Newtonian solvent. The dumbbells comprise two beads connected by a spring. The behavior of the dumbbell can be completely described by the vector $\mathbf{Q}$ that connects the location of one of the beads to the other. The governing equations are expressed solely in terms of $\mathbf{Q}$. To form the equations that describe the evolution of $\mathbf{Q}$, an elastic law for the spring must be proposed. We adopt the FENE relation for the spring force ${\mathbf{F}}^c$,

$${\mathbf{F}}^c=\frac{H\mathbf{Q}}{1-Q^2/Q_0^2},$$

(1)

where H is a spring constant and $Q_0$ is the maximum possible extension of the dumbbells [1]. Taking the limit $Q_0\to \infty$, the linear Hookean spring law is recovered, which microscopically reduces to the Oldroyd B model. In the limit at which the end-to-end vector reaches its maximum extension, the fluid acts analogous to a suspension of infinitely thin rigid rods. (1) is an empirical approximation to the physically derived inverse Langevin operator [25]. Other approximations could have been used, for instance, the Marko–Siggia force law [26], the Padé law [27], or improved closure approximation modifications such as the FENE-L model [28,29].

The equation of motion for an individual dumbbell can be found from the spring law. Assuming that the velocity field is homogeneous over the length scale of the polymer, then the evolution of $\mathbf{Q}$ is given by the Langevin equation [30]:

$$\mathbf{dQ}=\left(\mathit{L}·\mathbf{Q}-\frac{2}{\zeta }\left({\mathbf{F}}^c+{\mathbf{F}}^{MF}\right)\right)dt+\sqrt{\frac{4kT}{\zeta }}\mathbf{dW}.$$

(2)

Here, ${\mathbf{F}}^{MF}$ denotes the mean-field forces, k is the Boltzmann constant, T is the absolute temperature, and $\zeta$ is the drag constant. The tensor $\mathit{L}$ is the transpose of the velocity gradient given by $L_{ij}=\frac{\partial u_i}{\partial x_j}$. The Brownian term $\mathbf{dW}$ is a Wiener process whose components are independent Gaussian variables with mean zero and variance $dt$. The Langevin equation has an equivalent Fokker–Planck (FP) equation for the probability distribution $\psi (\mathbf{Q},t)$. The evolution equation for $\psi$ is given by the partial differential equation

$$\frac{\partial \psi }{\partial t}=\nabla ·\left(-\mathit{L}\mathbf{Q}\psi +\frac{2}{\zeta }\left({\mathbf{F}}^c+{\mathbf{F}}^{MF}\right)\psi \right)+\frac{2kT}{\zeta }{\nabla }^2\psi .$$

(3)

The velocity field is given by

$$\mathbf{u}=\left(\dot{\epsilon }x,-\dot{\epsilon }y\right),$$

(4)

where $\dot{\epsilon }$ is the elongation rate. We adopt the mean-field force term used in [12,31],

$${\mathbf{F}}^{\left(MF\right)}=-\frac{kT}{Q_0^4}f〈\stackrel{﹇}{\mathbf{QQ}}〉\mathbf{Q}.$$

(5)

The mean-field term is introduced to model weak anisotropic effects from concentration effects due to intra- and intermolecular interactions. The parameter f determines the strength of the mean-field term and is a positive function of the concentration. For infinitely dilute solutions, $f=0$. f increases with increasing concentration by $f={\left(c/c^{*}\right)}^{\frac{1}{3}}$, where $c^{*}$ is a reference concentration. $\stackrel{﹇}{\mathbf{QQ}}$ denotes the the deviatoric part of the symmetric decomposition of a tensor $\mathbf{QQ}$. Equation (5) was derived from studies of rod-like polymers [32]. At equilibrium, when $\mathbf{Q}\mathbf{Q}$ is isotropic, there is no preferential direction and the mean-field term becomes zero.

For finite $Q_0$, the system cannot be solved analytically as a closed constitutive equation cannot be formed from only the second moment, or radius of gyration tensor, $〈\mathbf{Q}\mathbf{Q}〉$. A closed system can be formed if the Hookean spring force is used; however, in elongational flows, one can expect the solution to be valid over a small range of elongation rates only due to the dumbbells becoming infinitely extended. The standard approach to overcome this closure problem is to replace the spring force model in (1) with the pre-averaged Peterlin closure approximation in which the nonlinear Wagner spring force term is approximated by a self-consistent average term [33]. This results in a spring force model of the form

$${\mathbf{F}}^c\approx \frac{H\mathbf{Q}}{1-〈Q^2〉/Q_0^2}.$$

(6)

This has the advantage that the nonlinear finite extensiblity term is now linear with respect to the microscopic variables.

3. Results

3.1. Mean-Field Force under Elongation

We now seek an exact solution to (3) for homogeneous planar elongational flow with elongation rate $\dot{\epsilon }$. This was first investigated in [13]. We extend their analysis to investigate the non-Hookean response. We will use the convention that $\mathit{L}$ is the velocity gradient tensor, i.e., $L_{ij}=\frac{\partial u_i}{\partial x_j}$. For purely planar elongational flow, $\mathit{L}$ is given by $L_{ij}=\dot{\epsilon }{\delta }_{i1}{\delta }_{j1}-\dot{\epsilon }{\delta }_{i2}{\delta }_{j2}$, where ${\delta }_{ij}$ is the Kronecker delta function. To find a steady-state solution, we adopt a change of variable approach by writing $\psi =\chi ·\varphi$, where $\varphi =e^{\frac{\zeta }{4kT}\mathbf{L}:\mathbf{Q}\mathbf{Q}}$ and: denotes full index contraction. This form of solution is motivated by solutions of the FP equation found for elongational flows in [30]. This removes the advective term and the FP equation reduces to

$$\nabla ·\left(\frac{2}{\zeta }\frac{H\mathbf{Q}\chi }{\left(1-〈Q^2〉/Q_0^2\right)}-\frac{2kT}{\zeta Q_0^4}f〈\stackrel{﹇}{\mathbf{QQ}}〉\mathbf{Q}\chi \right)+\frac{2kT}{\zeta }{\nabla }^2\chi =0.$$

(7)

As we are concerned with pure elongational flow, we can assume a priori that the shear term ${〈\mathbf{QQ}〉}_{12}=0$. Thus, the tensor product term in (5) can be written as

$$〈\stackrel{﹇}{\mathbf{QQ}}〉{}_{ip}{\mathbf{Q}}_p=〈\stackrel{﹇}{\mathbf{QQ}}〉{}_{\left[ii\right]}{\mathbf{Q}}_{\left[i\right]},$$

(8)

where $\left[\right]$ implies no summation over the elements. It follows that at the microscopic level the spring force is now a scalar multiplied by the end-to-end vector. This additional force acts analogous to a Hookean spring with non-constant H, where the spring constant is a function of the macroscopic variables. Thus we can write a solution for the probability distribution $\psi$ as

$$\begin{array}{c}\psi (r,\theta )=Jexp\left\{-\frac{1}{2}W{r}^2\left(1-\epsilon \left({cos}^2\left(\theta \right)-{sin}^2\left(\theta \right)\right)\right)\right\}exp\left\{\frac{H\lambda \dot{\epsilon }{r}^2}{kT}({cos}^2\left(\theta \right)-{sin}^2\left(\theta \right))\right\}.\end{array}$$

(9)

To simplify the notation, we replace $Q_r$, $Q_{\theta }$, $Q_x$, $Q_y$ by $r,\theta ,x$, y and introduce the characteristic spring time, $\lambda =\zeta /4H$. The variables $W,J$ and $\epsilon$ are given by

$$\begin{array}{c}W=\frac{H}{kT\left(1-〈r^2〉/Q_0^2\right)},J=\frac{H\sqrt{\left\{1-\epsilon -2\lambda \dot{\epsilon }\left(1-〈r^2〉/Q_0^2\right)\right\}\left\{1+\epsilon +2\lambda \dot{\epsilon }\left(1-〈r^2〉/Q_0^2\right)\right\}}}{2\pi kT\left(1-〈r^2〉/Q_0^2\right)},\\ \epsilon =\frac{kT}{H{Q}_0^4}f〈\stackrel{﹇}{\mathbf{QQ}}〉{}_{11}\left(1-〈Q{}^2〉/Q_0^2\right).\end{array}$$

(10)

Physically $\epsilon$ is a measure of the anisotropy that arises solely from mean-field effects and acts like an additional elongation rate. $〈r^2〉$ and $\epsilon$ are not free variables as they have a hidden dependency on $\psi$. They can be determined by the averaging conditions

$$\begin{array}{c}〈r^2〉={\int }_0^{2\pi }{\int }_0^{\infty }{r}^2\psi (r,\theta ;〈r^2〉,\epsilon )rdrd\theta ,\end{array}$$

(11)

$$\begin{array}{c}\frac{2H{Q}_0^4\epsilon }{kTf\left(1-〈r^2〉/Q_0^2\right)}={\int }_0^{2\pi }{\int }_0^{\infty }{r}^2\left({cos}^2\left(\theta \right)-{sin}^2\left(\theta \right)\right)\psi (r,\theta ;〈r^2〉,\epsilon )rdrd\theta .\end{array}$$

(12)

Integration of (11) and (12) gives the transcendental equations

$$\begin{array}{c}〈r^2〉=\frac{2kT\left(1-〈r^2〉/Q_0^2\right)}{H(1-\epsilon -2\lambda \dot{\epsilon }\left(1-〈r^2〉/Q_0^2\right))(1+\epsilon +2\lambda \dot{\epsilon }\left(1-〈r^2〉/Q_0^2\right))},\end{array}$$

(13)

$$\begin{array}{c}\frac{H{Q}_0^4\epsilon }{kTf\left(1-〈r^2〉/Q_0^2\right)}=\frac{kT\left(\epsilon +2\lambda \dot{\epsilon }\left(1-〈r^2〉/Q_0^2\right)\right)\left(1-〈r^2〉/Q_0^2\right)}{H(1-\epsilon -2\lambda \dot{\epsilon }\left(1-〈r^2〉/Q_0^2\right))(1+\epsilon +2\lambda \dot{\epsilon }\left(1-〈r^2〉/Q_0^2\right))}.\end{array}$$

(14)

Adopting the scalings ${\mathbf{Q}}^{*}=\frac{\mathbf{Q}}{\sqrt{KT/H}}$, $t^{*}=\frac{t}{\lambda },{\dot{\epsilon }}^{*}=\lambda \dot{\epsilon }$, ${\psi }^{*}=\frac{kT}{H}\psi$ and $f^{*}=\frac{k^2{T}^2}{H^2{Q}_0^4}f$, we obtain

$$\begin{array}{c}{\widehat{\mathbf{F}}}^c=\frac{1}{2}\frac{{\mathbf{Q}}^{*}}{1-〈Q^2〉/b},{\widehat{\mathbf{F}}}^{MF}=-\frac{1}{2}{f}^{*}〈\stackrel{﹇}{\mathbf{Q}*\mathbf{Q}*}〉{\mathbf{Q}}^{*},\\ \frac{\partial {\psi }^{*}}{\partial t^{*}}=\nabla ·\left(-{\mathit{L}}^{*}{\mathbf{Q}}^{*}{\psi }^{*}+\left({\widehat{\mathbf{F}}}^c+{\widehat{\mathbf{F}}}^{MF}\right){\psi }^{*}\right)+\frac{1}{2}{\nabla }^2{\psi }^{*}.\end{array}$$

(15)

The parameter, $b=H{Q}_0^2/kT$, is the ratio of the flexible spring time $\frac{\zeta }{4H}$ to the characteristic rigid time $\frac{\zeta L^2}{12kT}$. Physically b represents the number of monomers in the polymer [34]. Dropping the ${}^{*}$ for convenience, we recover a system of nonlinear algebraic equations:

$$\begin{array}{c}〈r{}^2〉=\frac{2\left(1-〈r{}^2〉/b\right)}{(1-\epsilon -2\dot{\epsilon }\left(1-〈r{}^2〉/b\right))(1+\epsilon +2\dot{\epsilon }\left(1-〈r{}^2〉/b\right))},\end{array}$$

(16)

$$\begin{array}{c}\epsilon =\frac{f\left(\epsilon +2\dot{\epsilon }\left(1-〈r{}^2〉/b\right)\right){\left(1-〈r{}^2〉/b\right)}^2}{(1-\epsilon -2\dot{\epsilon }\left(1-〈r{}^2〉/b\right))(1+\epsilon +2\dot{\epsilon }\left(1-〈r{}^2〉/b\right))}.\end{array}$$

(17)

3.2. Mean-Field Force under Shear

The effects of concentration dependent shear-induced anisotropy were investigated in [12] by means of solving the constitutive equation for the second moment. However, at no extra computational cost, we can obtain $\psi$ which permits the computation of any macroscopic variable, thus enabling direct comparison of the predictions of the model with molecular dynamics simulations.

In contrast to the case for elongational flow, it is clear that the term ${〈\mathbf{QQ}〉}_{12}$ will be non-zero and thus a different approach must be taken. For Hookean dumbbells under shear flow with no mean-field dependence, the velocity gradient tensor $\mathit{L}$ is given by $L_{ij}=\dot{\gamma }{\delta }_{i1}{\delta }_{j2}$, where ${\delta }_{ij}$ is the Kronecker delta function, and the FP equation has solution [30] takes the form

$${\psi }_{SR}\propto exp\left\{-\frac{H}{2kT}\frac{\left(x^2-2\lambda \dot{\gamma }xy+(1+2{\lambda }^2{\dot{\gamma }}^2)y^2\right)}{1+{\lambda }^2{\dot{\gamma }}^2}\right\}.$$

(18)

For the analogous case, which includes a mean-field force under shear, we anticipate that the solution is of a similar form as (18), and seek a trial solution of the form

$$exp\left(-\frac{H(A{x}^2+Bxy+C{y}^2)}{2kT(1-〈Q^2〉/Q_0^2)}\right).$$

(19)

where $A,B,C$ are constants to be determined. We need not restrict ourselves to pure shear flow but could consider a flow field that exhibits both elongational and shear phenomena. We take a general homogeneous velocity gradient field for which the velocity components are given by $u=\dot{\epsilon }x+{\dot{\gamma }}_{1}y$ and $v=-\dot{\epsilon }y+{\dot{\gamma }}_{2}x$. Substitution of this trial solution into the FP equation leads to

$$\nabla ·\left[-\dot{\epsilon }\mathit{q}\psi -\left(\begin{array}{c}{\dot{\gamma }}_1\\ {\dot{\gamma }}_2\end{array}\right)\circ {\mathit{q}}^{\prime }\psi +\frac{2H\mathbf{Q}\psi }{\zeta (1-〈Q^2〉/Q_0^2)}+\frac{2kT}{\zeta }\nabla {\psi }_{}\right]=0$$

(20)

where $\mathit{q}=(x,-y)$ and ${\mathit{q}}^{\prime }=(y,x)$. The symbol ∘ denotes array-wise multiplication. We find

$$\begin{array}{c}A=\frac{(1-2{\dot{\gamma }}_1{\dot{\gamma }}_2{\lambda }^2+2{\lambda }^2{\dot{\gamma }}_2^2-2\lambda \epsilon )}{(1+{\dot{\gamma }}_1^2{\lambda }^2-2{\lambda }^2{\dot{\gamma }}_1{\dot{\gamma }}_2+{\dot{\gamma }}_2^2{\lambda }^2)},B=\frac{-2\lambda \left(({\dot{\gamma }}_1+{\dot{\gamma }}_2)+2\epsilon \left({\dot{\gamma }}_2-{\dot{\gamma }}_1\right)\lambda \right)}{(1+{\dot{\gamma }}_1^2{\lambda }^2-2{\lambda }^2{\dot{\gamma }}_1{\dot{\gamma }}_2+{\dot{\gamma }}_2^2{\lambda }^2)},\\ C=\frac{(1-2{\lambda }^2{\dot{\gamma }}_1{\dot{\gamma }}_2+2{\lambda }^2{\dot{\gamma }}_1^2+2\lambda \epsilon )}{(1+{\dot{\gamma }}_1^2{\lambda }^2-2{\lambda }^2{\dot{\gamma }}_1{\dot{\gamma }}_2+{\dot{\gamma }}_2^2{\lambda }^2)}.\end{array}$$

(21)

The mean-field term is decomposed into a shear-like component and an elongation-like component: ${\mathbf{F}}^{\left(MF\right)}=-\frac{kT}{2Q_0^4}f\left(〈x^2-y^2〉\right)\mathit{q}-\frac{kT}{Q_0^4}f〈xy〉{\mathit{q}}^{\prime }$, by means of two basis vectors, $\mathit{q}$ and ${\mathit{q}}^{\prime }$. Therefore,

$$\nabla ·\left[-\frac{kT}{\zeta Q_0^4}f\left(〈x^2-y^2〉\right)\mathit{q}\psi -\left(\begin{array}{c}\dot{\gamma }+\frac{2kT}{\zeta Q_0^4}f〈xy〉\\ \frac{2kT}{\zeta Q_0^4}f〈xy〉\end{array}\right)\circ {\mathit{q}}^{\prime }\psi +\frac{2H\mathbf{Q}\psi }{\zeta (1-〈Q^2〉/Q_0^2)}+\frac{2kT}{\zeta }\nabla {\psi }_{}\right]=0,$$

(22)

which is analogous to the solution for $f=0$ with modified shear rate ${\dot{\gamma }}_1=\dot{\gamma }+{\dot{\gamma }}_2,{\dot{\gamma }}_2=\frac{2kTf}{\zeta Q_0^4}〈xy〉$ and modified elongation rate $\dot{\epsilon }=\frac{kTf}{\zeta Q_0^4}〈x^2-y^2〉$. From (21) the parameters $A,B,C$ are

$$\begin{array}{c}A=\frac{1-\dot{\gamma }\lambda \mathsf{\Gamma }{\left(1-〈r^2〉/Q_0^2\right)}^2-\mathsf{\Theta }\left(1-〈r^2〉/Q_0^2\right)}{1+{\lambda }^2{\dot{\gamma }}^2{\left(1-〈r^2〉/Q_0^2\right)}^2},\\ B=\frac{-2\dot{\gamma }\lambda \left(1-〈r^2〉/Q_0^2\right)-2\mathsf{\Gamma }\left(1-〈r^2〉/Q_0^2\right)+2\dot{\gamma }\lambda \mathsf{\Theta }{\left(1-〈r^2〉/Q_0^2\right)}^2)}{1+{\lambda }^2{\dot{\gamma }}^2{\left(1-〈r^2〉/Q_0^2\right)}^2},\\ C=\frac{1+2{\dot{\gamma }}^2{\lambda }^2\left(1-〈r^2〉/Q_0^2\right)+\dot{\gamma }\lambda \mathsf{\Gamma }{\left(1-〈r^2〉/Q_0^2\right)}^2+\mathsf{\Theta }\left(1-〈r^2〉/Q_0^2\right)}{1+{\lambda }^2{\dot{\gamma }}^2{\left(1-〈r^2〉/Q_0^2\right)}^2},\end{array}$$

(23)

where $\mathsf{\Theta }=\frac{kTf}{H{Q}_0^4}〈x^2-y^2〉/2$ and $\mathsf{\Gamma }=\frac{kTf}{H{Q}_0^4}〈xy〉$. In the limit as f tends to zero (i.e., $\mathsf{\Theta },\mathsf{\Gamma }\to 0$) and $Q_0\to \infty$ the standard Hookean solution (18) is recovered. The normalization constant J can be expressed as $J=\frac{H\sqrt{4AC-B^2}}{4\pi kT\left(1-〈r^2〉/Q_0^2\right)}$. The consistency conditions for the system naturally give rise to an equivalent set of equations as derived by [12] who used a change of basis method. The two conditions (11) and (12) still hold, however, as no ansatz is used regarding $A_{12}$, the extra degree of freedom is accounted for by an additional consistency condition. We can write

$$\begin{array}{c}〈r^2〉=J\int r^2\psi (r,\theta ;〈r^2〉,\mathsf{\Gamma },\mathsf{\Theta })rdrd\theta ,\end{array}$$

(24)

$$\begin{array}{c}\frac{H{Q}_0^4}{kTf}\mathsf{\Gamma }=J{\int }_0^{2\pi }{\int }_0^{\infty }{r}^{2}sin\left(\theta \right)cos\left(\theta \right)\psi (r,\theta ;〈r^2〉,\mathsf{\Gamma },\mathsf{\Theta })rdrd\theta ,\end{array}$$

(25)

$$\begin{array}{c}\frac{H{Q}_0^4}{kTf}\mathsf{\Theta }=\frac{1}{2}J{\int }_0^{2\pi }{\int }_0^{\infty }{r}^2\left({cos}^2\left(\theta \right)-{sin}^2\left(\theta \right)\right)\psi (r,\theta ;〈r^2〉,\mathsf{\Gamma },\mathsf{\Theta })rdrd\theta .\end{array}$$

(26)

Use of the aforementioned scalings, with the shear rate scaling analogous to the elongation rate (i.e., ${\dot{\gamma }}^{*}=\dot{\gamma }\lambda$ ), leads to the system of equations (with the star notation dropped),

$$\begin{array}{c}〈r^2〉=\frac{2\left(1+{\dot{\gamma }}^2{\left(1-〈r^2〉/b\right)}^2\right)\left(1-〈r^2〉/b\right)}{\left(1-{\left(1-〈r^2〉/b\right)}^2\left({\mathsf{\Gamma }}^2+{\mathsf{\Theta }}^2+2\dot{\gamma }\mathsf{\Gamma }\right)\right)},\end{array}$$

(27)

$$\begin{array}{c}\mathsf{\Gamma }=\frac{f\left(\dot{\gamma }+\mathsf{\Gamma }-\dot{\gamma }\mathsf{\Theta }\left(1-〈r^2〉/b\right)\right){\left(1-〈r^2〉/b\right)}^2}{\left(1-{\left(1-〈r^2〉/b\right)}^2\left({\mathsf{\Gamma }}^2+{\mathsf{\Theta }}^2+2\dot{\gamma }\mathsf{\Gamma }\right)\right)},\end{array}$$

(28)

$$\begin{array}{c}\mathsf{\Theta }=\frac{f\left({\dot{\gamma }}^2\left(1-〈r^2〉/b\right)+\mathsf{\Theta }+\dot{\gamma }\mathsf{\Gamma }\left(1-〈r^2〉/b\right)\right){\left(1-〈r^2〉/b\right)}^2}{\left(1-{\left(1-〈r^2〉/b\right)}^2\left({\mathsf{\Gamma }}^2+{\mathsf{\Theta }}^2+2\dot{\gamma }\mathsf{\Gamma }\right)\right)},\end{array}$$

(29)

which are again solved by the Newton–Raphson method. It can be seen that for elongational flow, increasing f leads to greater stretching of the dumbbell. For shear flows, increasing f leads to greater orientation of the dumbbell along the flow. We will now compare these results to those of the full nonlinear FENE model. The expansion for small $\dot{\gamma }$ was previously obtained in [12], however the effect for large shear was not addressed. We will consider this effect by using the method of dominant balance [35]. This results in the expressions

$$\begin{array}{c}〈r^2〉\sim b-2^{-\frac{1}{3}}{b}^{\frac{4}{3}}{\dot{\gamma }}^{-\frac{2}{3}}+\frac{2^{\frac{1}{3}}{b}^{\frac{2}{3}}\left(2+b+b^{2}f\right)}{6}{\dot{\gamma }}^{-\frac{4}{3}}+\mathcal{O}\left({\dot{\gamma }}^{-\frac{6}{3}}\right)\\ \mathsf{\Theta }\sim \frac{bf}{2}-\frac{2^{\frac{2}{3}}{b}^{\frac{1}{3}}f\left(2+b\right)}{4}{\dot{\gamma }}^{-\frac{2}{3}}+\mathcal{O}\left({\dot{\gamma }}^{-\frac{4}{3}}\right),\\ \mathsf{\Gamma }\sim 2^{-\frac{2}{3}}{b}^{\frac{2}{3}}f{\dot{\gamma }}^{-\frac{1}{3}}-\frac{f}{12}\left(b^{2}f+4b+8\right){\dot{\gamma }}^{-1}+\mathcal{O}\left({\dot{\gamma }}^{-\frac{4}{3}}\right).\end{array}$$

These series exhibit much better convergence than those for elongational flow. We see that the effects of concentration increase the radius of gyration. These effects decay as ${\dot{\gamma }}^{-\frac{4}{3}}$, which is a slower decay than for the elongational case, where the decay is proportional to ${\dot{\epsilon }}^{-2}$. An interesting feature of shear flows is that the degree of shear changes the orientation of the dumbbells. The orientation angle is given by [36]

$$\chi =\frac{\pi }{4}-\frac{1}{2}{tan}^{-1}\left(\mathsf{\Gamma }{\mathsf{\Theta }}^{-1}\right)\sim \frac{\pi }{4}-2^{-\frac{2}{3}}{b}^{\frac{2}{3}}{\dot{\gamma }}^{-\frac{1}{3}}+\frac{1}{12}\left(bf-2\right){\dot{\gamma }}^{-1}.$$

(30)

As $\dot{\gamma }\to \infty$, the angle of orientation is $\frac{\pi }{4}$ as expected and the dominant effect is the finite extension force. We see that the mean-field term which decays as ${\dot{\gamma }}^{-1}$ has opposite sign to the finite extension effect and acts so as to resist the orientation of the dumbbell.

3.3. FENE Solution

An iterative finite element scheme was used to solve the FP equation with the inclusion of the fully nonlinear FENE force term. Although the Peterlin closure approximation removes the unwanted infinite elongation that occurs with Hookean spring models, it only imposes a constraint on the averaged value of the end-to-end vector rather than an upper-bound, thus dumbbells can exist which have length greater than $Q_0$ [37]. Furthermore, the effect of the Peterlin closure approximation is non-local. The distribution for the FENE-P model is Gaussian whereas the FENE model cannot have a chain length greater than the maximum extension that gives compact support. This difference in the FENE-P model has previously been shown to effect the transient response as well [14]. In this study, we limit our investigation to steady state flows. We will use the notation that $\mathsf{\Omega }$ is the computational domain and $\partial \mathsf{\Omega }$ is the bounding region. Under purely elongational flow, the FP equation is solved for radius $0<r<\sqrt{b}$, which arises from the FENE model restricting the radius of gyration. To reduce the computational domain, we consider the upper quadrant and impose symmetry boundary conditions along the $x-$ and $y-$axes. The Dirichlet condition, $\psi =0$, is applied on the boundary $r=\sqrt{b}$, which enforces zero probability of the dumbbell being fully extended. Strong symmetry does not exist for the pure shear case and the FP equation must be solved over the entire circular domain. As such a system has purely homogeneous boundary conditions, the trivial zero solution would be recovered. To overcome this we add a point constraint on the probability at the origin where the constant is determined by requiring conservation of probability.

Simulations from FENE and FENE-P schemes for elongational flow are shown in Figure 1 and Figure 2 for $f=0$, 0.5 and 0.75. The analytical solution for $f=0$ is given by $N\int r^2{(1-r^2/b)}^{b/2}{e}^{\dot{\epsilon }{r}^{2}cos\left(2\theta \right)}d\mathsf{\Omega }$ where N is the appropriate normalization constant [7]. This is in close agreement with the results from the finite element solution (Figure 2a). Plots of the radius of gyration calculated from the FENE and FENE-P models for $f=0$, 0.5, and 0.75 for the case of shear flow are shown in Figure 3. It was found that for large elongation rates, use of the Peterlin closure approximation (Figure 2) gave similar results to the FENE model for the dependence of the radius of gyration on the non-dimensionalized elongation rate. However, the Peterlin closure approximation systematically overestimated the radius of gyration. For the case $f=0$, the overestimation can be explained as the finite extension is imposed only on the average extension. We also found that the greatest error occurs for moderately small elongation rates. The magnitude of the error is increased and the maximum peak error occurs at a lower elongation rate for increasing values of f. This error is due to the relatively large gradients of the radius of gyration with respect to the elongation rate. These large gradients are due to coil stretch transition and are eventually smoothed when the finite extensibility constraint starts to dominate.For shear flow with $f=0$, the error increases monotonically with shear rate. However, with the introduction of the mean-field term the peak error occurs at a finite shear rate. The effect of increasing f is similar, but less pronounced than the effect seen with elongational flows, where the magnitude of the error increases and occurs at a lower shear rate.

Using the previous scaling, while dropping the ${}^{*}$, the steady-state FP equation with the inclusion of a mean-field term for elongational flow is given by

$$\nabla ·\left[-\dot{\epsilon }(x{\mathit{e}}_{\mathit{x}}-y{\mathit{e}}_{\mathit{y}})\psi +\frac{1}{2}\left(\frac{\mathit{x}}{1-r^2/b}-f〈\stackrel{﹇}{\mathbf{QQ}}〉\mathit{x}\right)\psi \right]+\frac{1}{2}{\nabla }^2\psi =0,$$

(31)

where ${\mathit{e}}_{\mathit{x}}$ and ${\mathit{e}}_{\mathit{y}}$ are the unit vectors in the x- and y-directions.

Equation (31) does not account for the following two constraints. Firstly, that ${\int }_{\mathsf{\Omega }}\psi dV=1$, and second, upon scaling of $\psi$, the condition of $〈\stackrel{﹇}{\mathbf{QQ}}〉=\int \stackrel{﹇}{\mathbf{QQ}}\psi dV$ which introduces nonlinearity in $\psi$. To approach these conditions a fixed point finite element iterative scheme is used to reach a default tolerance:

$$\begin{array}{c}\mathcal{L}(\psi ,\mathsf{\Phi };\dot{\epsilon },〈\stackrel{﹇}{\mathbf{QQ}}〉{}_n)=0,\end{array}$$

(32)

$$\begin{array}{c}{\psi }^{*}=\frac{\psi }{{\int }_{\mathsf{\Omega }}\psi dS},\end{array}$$

(33)

$$\begin{array}{c}〈\stackrel{﹇}{\mathbf{QQ}}〉{}_{n+1}={\int }_{\mathsf{\Omega }}\stackrel{﹇}{\mathbf{QQ}}{\psi }^{*}d\mathbf{Q},\end{array}$$

(34)

where $\mathcal{L}\left(\psi ,\mathsf{\Phi };\dot{\epsilon },〈\stackrel{﹇}{\mathbf{QQ}}〉\right)$ denotes the weak form of (31). Equation (31) can be similarly formed and solved for the case of pure shear flow which in weak form is

$$\begin{array}{c}\hfill {\int }_{\mathsf{\Omega }}\nabla ·\left(-\mathit{L}\mathit{x}\psi -\frac{1}{2}\left(\frac{\mathit{x}}{1-r^2/b}-f〈\stackrel{﹇}{\mathbf{QQ}}〉\mathit{x}\right)\psi \right)\mathsf{\Phi }-\frac{1}{2}\nabla \psi ·\nabla \mathsf{\Phi }dS+\frac{1}{2}{\int }_{\partial \mathsf{\Omega }}\mathsf{\Phi }\nabla \psi ·\mathit{dr}=0,\end{array}$$

where $\mathsf{\Phi }$ gives the shape function. The velocity gradient tensor, $\mathit{L}$, takes a different form for shear and elongational flows. The Galerkin finite element method was used, where $\mathsf{\Phi }$ comprised linear piecewise polynomials. The stiffness matrix was formed from the first term which is integrated using a fourth order quadrature method over a triangular mesh with 40,000 nodes. The mesh density was increased near $r=b$, as well as along $y=0$ for the elongational case, to capture the larger gradients in $\psi$. The advective term is the only term which changes and is initially set to be 0. ${〈\mathbf{QQ}〉}_{12}$ is found through the iterative scheme (32)–(34). The probability density function $\psi$ is plotted for the FENE solution for the case of pure elongational flow in Figure 4. The effects of the additional concentration are similar to those obtained using the Peterlin closure approximation where the concentration term acts analogously to the elongation term and gives rise to increased probability of larger extensions of the dumbbells. The effect is similar, but smaller, for the case of pure shear flow as shown in Figure 5.

3.4. Comparison with Molecular Dynamics Simulations

The predictions of the FENE-P model were compared with stochastic data [23]. The FENE-P model predicts that with increasing f, the mean-field effect will increase the mean extension of the dumbbell. It can be seen from (7) and (8) that the spring stiffness constant reduces for increasing f. This feature is opposite to that which is reported in the literature. One possible explanation is that in the model used in [12] the equilibrium radius of gyration is independent of concentration. However, the authors of [23] report that the equilibrium radius of gyration can vary with concentration. We model the equilibrium mean-field effects by expressing the parameter b as a function of the concentration. We model b as a quadratic in $c/c^{*}$, i.e.,

$$b=b_e\left(1+b_1\left(\frac{c}{c^{*}}\right)+b_2{\left(\frac{c}{c^{*}}\right)}^2\right),$$

where the parameters $b_1,b_2$ are found by fitting to the profile of the equilibrium results for ${10}^{-3}<\frac{c}{c^{*}}<10$ given in [23]. Stoltz et al. [23] determined the extension as a function of both the relaxation time and the concentration dependent relaxation time (CDRT). The CDRT is found by modeling the experimental procedure in [18], in which the polymer was stretched under elongational flow and then allowed to relax back to equilibrium. The CDRT was found by fitting an exponential function to the transient profile of the radius of gyration. Stoltz et al. [23] record the affect of concentration on the mean distance between the end points only in streamwise direction and not the radius of gyration. This can be found using probability density function (9) to be

$$〈X〉={\left\{\frac{2\left(1-〈Q^2〉/b\right)}{\pi \left[1-\epsilon +2\dot{\epsilon }\left(1-〈Q^2〉/b\right)\right]}\right\}}^{\frac{1}{2}}.$$

(35)

We chose $b_e=100$ for our simulations as the results were found to be largely independent of $b_e$ provided it was sufficiently large. Figure 6 shows plots of $〈X〉$ normalized by its value for large extension against the elongation rate $\dot{\epsilon }$ and against the critical elongation rate ${\dot{\epsilon }}_c$. The model now predicts that the end-to-end distance is reduced for increasing concentration when plotted against ${\dot{\epsilon }}_c$. For comparison, $〈X〉$ is plotted in Figure 7, for the case where the mean-field force term was set to zero. If we look at the normalized values of $\dot{\epsilon }$$〈X〉$ corresponding to $\dot{\epsilon }=1$ we observe a small overshoot in for the case where the mean-field term is included (Figure 6a). Such an overshoot does not occur in the simulations where this term has been neglected (Figure 7). A similar phenomena was reported in [23], although in their study, the size of the overshoot was more pronounced.

4. Discussion

The principle result of this article was the derivation of solutions to the Fokker–Planck equation incorporating a FENE-P dumbbell model with the addition of a mean-field term for both shear and elongational flows, including the effects of strong shear and elongations close to the maximum permissible. We found that the results corresponding to elongational flow were amenable to analytical analysis. For a linear spring force, the distribution can be found exactly and for dumbbells near full extension, an exact solution was found that could be expressed as the root of a quartic equation. The solution for a constant shear flow was found by directly solving the Fokker–Planck equation. We further investigated the effect that the closure assumption on the spring force has on the model. The probability density function was found numerically using a finite element scheme. We found that the closure assumption in conjunction with the mean-field term led to an over-estimation of the extension of the dumbbell at moderate values of the elongation and shear rate. In the latter part of the work, we matched the prediction of the model to molecular dynamics simulations. It was found that the model developed by [12] fails to predict the reduced stretching that occurs with increasing concentration. We thus proposed a novel extension by introducing a concentration dependence on the parameter b, which describes the ratio of the flexible spring time to the characteristic rigid time, to match the results given in [23], which were based on stochastic data. This modified model correctly predicted reduced stretching with increased concentration; however, the overshoot in relation to the infinitely dilute case in the aligned end-to-end vector was underpredicted.

Throughout this article we have only considered mean-field effects under steady shear and steady elongation conditions. This leaves the areas of transient shear and elongation open to further study.