On the similarities between the simplified Phan-Thien–Tanner model and the finitely extensible nonlinear elastic dumbbell (Peterlin closure) model in simple and complex flows

Of the many important contributions to rheology and non-Newtonian fluid mechanics by Professor R. B. Bird that have passed the test of time, we choose in this commemorative special issue paper to highlight the constitutive model for dilute solutions of polymers based on kinetic theory, now known (following Bird's suggestion) by the widespread acronym of the “FENE-P” model.^1–4 In particular, we focus on the connections to the seemingly unsimilar model due to Phan-Thien and Tanner,⁵ which was originally developed for polymeric melts and concentrated solutions using network theory.^6,7 Bird's interest in molecular theory (or kinetic theory as he preferred to call it) seems to have been motivated by a desire to “provide a clue as to which forms of constitutive equations may prove to be the most useful.”⁸ [In a fascinating oral history,⁹ Bird recounts that on first hearing a talk on the kinetic theory of rigid dumbbells “at that time I thought, ‘What a lot of nonsense. Why would anybody be interested in rigid dumbbells?’ (A point of view that I have changed very much since).”] Starting in the early 1970s with a “survey” paper¹⁰ written with a student (Hal Warner) and a post doc (D. Colin Evans—a former Ph.D. student of J. G. Oldroyd), and initially focusing on rigid macromolecules (idealized as rigid dumbbells), Bird and collaborators^11–13 made significant strides in developing constitutive models derived from kinetic theory culminating in the classic two-volume monograph “Dynamics of Polymeric Liquids”^2,3 (referred to as “DPL” hereafter) with the second volume dedicated specifically to this topic. Bird's continued his contributions to the subject between the first and second editions of DPL, especially with regard to kinetic theory of finitely extensible (“FENE”) dumbbells, such that the second edition must have required significant revision! Using kinetic theory to understand rheological behavior, and in particular, the “FENE-P” dumbbell model, which we focus on here, held a long fascination for Professor Bird.^14–16 The benefit of the kinetic theory approach, particularly from a pedagogical perspective, was something which obviously appealed.¹⁷ 

In this special issue paper, we compare the response of this classical FENE-P dumbbell model^1–4 derived for dilute polymeric solutions with the simplified form of the linear Phan-Thien–Tanner fluid model⁵ derived from network theory (for polymeric melts and concentrated solutions). Although the potential similarity between these models was highlighted even in the original paper of Bird et al.¹ [“Their equation is similar to our Eqs. (57) and (60) and exhibits the same kinds of nonlinearities in the stresses.”] and has been discussed in passing in a number of papers in the intervening 40 years,^18–22 here we probe the two models across a significantly broader range of conditions and flows than has been done hitherto. In so doing, we highlight, despite their completely different microstructural origin, the very similar behavior of the two models. Before describing the results of our investigations, we first provide an historical overview of the development of the two models. We then discuss the response of the two models in flows of increasing complexity from simple homogeneous in space and time, to homogeneous in space but unsteady in time, to unsteady in both space and time. Finally, we compare the response of the two models in a range of complex geometries, which although restricted to Eulerian steady flow, are unsteady in a Lagrangian sense.

The FENE-P dumbbell model was derived considering the dynamics of a basic dumbbell molecule formed by two spheres connected by a spring with end-to-end vector $Q$ (a stochastic quantity) subject to the influence of Brownian thermal motion, viscous drag, and an elastic relaxation by a non-linear Warner spring,²³ resulting in an equation for the evolution of the first moment $⟨QQ⟩$ and another to obtain stress.¹ Although those equations contained all that was necessary to calculate the flow and many elements of the molecular configuration (stretch, orientation), it was common at the time (in the 1980s) to transform the set into a stress formulation equation. In the case of the FENE-P dumbbell model, the first published reference was, in fact, by Bird,⁴ but the model was described more fully by Bird et al.¹ who gave

where the maximum molecular stretch $Q0$ (⁠$Q02≡⟨Q⋅Q⟩max$⁠) was normalized as $b=HQ02/kT$ (⁠$H$ is the Hookean spring constant, $k$ is the Boltzmann constant, $T$ is the absolute temperature, and $n$ is the number of molecules per unit volume). Here, $τ$ represents the polymeric contribution to stress, $I$ is the identity tensor, $γ̇$ is the rate of the deformation tensor (⁠$γ̇=∇u+∇uT$⁠), $λ$ is the Rouse relaxation time, and $(τ)∇$ represents the so-called upper-convected derivative of the stress tensor due to Oldroyd,²⁴ which is defined as

and $D/Dt$ is the usual material (or substantial) derivative [$Dx/Dt=∂x/∂t+(u.∇)x$]. At equilibrium, the polymer stress vanishes, the trace $τkk$ is, therefore, zero $τkk=τxx+τyy+τzz=0$⁠, so $Zeq=(b+3)/b$⁠, which we will denote simply by $a≡Zeq$⁠. Interestingly, in this paper, Bird et al.¹ stated that “It is appropriate to call Eq. (10) [i.e., Eqs. (1a) and (1b) here] Tanner's equation, since Tanner exactly obtained this same result, albeit by a different route.”²⁵ 

In another major contribution of Professor Bird, namely, the two volume set of Dynamics of Polymeric Liquids^2,3 (together with Armstrong, Hassager, and, for volume two, Curtiss), which constitutes a comprehensive treatise of non-Newtonian fluid mechanics and kinetic theory modeling of macromolecular fluids, Bird proposed an improvement to the original FENE-P dumbbell model such that the zero-shear-rate viscosity of the model coincided with that derived from the FENE theory²³ without the Peterlin²⁶ approximation, giving $ηp=nkTλb/(b+5)$ [in the 1980 paper, it was $ηp=nkTλb/(b+3)$⁠; note that $ηp$ is the polymeric viscosity $ηp=η(γ̇→0)−ηs$⁠, where $ηs$ is the solvent viscosity]. Such an improvement was achieved via the addition of an isotropic term to the standard Peterlin approximation, namely, $⟨QQ1−Q2/Q02⟩≅⟨QQ⟩1−⟨Q2⟩/Q02+eQ02I$⁠, and from calculation at equilibrium of the trace of this expression taken as an equality, the constant $e$ was found to be $e=2/(b(b+2))$⁠; see p. 89 of DPL2.² The original FENE-P dumbbell model of Bird et al.¹ in the stress formulation, after inserting the expressions for $aηp$ to replace $nkTλ$⁠, can be transformed to the improved version simply by setting $b→b+2$⁠. (Note that this rule does not apply to the definitions of $ηp$⁠.)

Nowadays, it is much more common to see the constitutive equations written in terms of the conformation formulation,^27–30 which ultimately are the basic equations when the principles of mechanics are applied in the context of kinetic theory, followed by phase averaging in a phase space and introducing approximations to arrive at a closed form manageable set of differential equations. We prefer to work with a conformation tensor that reduces to the identity tensor at equilibrium (a situation of no flow or of zero deformation rates), namely, $A=⟨QQ⟩/13⟨Q⋅Q⟩eq$⁠, and in this case, the previous equation for the FENE-P dumbbell model becomes a pair of equations for the evolution of $A$ and for a definition of stress $τ$ from the conformation,

and $L2=b+3$ in the original model¹ or $L2=b+5$ in the improved version of Bird et al.^2,3$L2$ has the physical meaning representing the square of the ratio between the fully extended dumbbell and its size at equilibrium [$L2=Q02/Req2$ where $Req2=13⟨Q⋅Q⟩eq$⁠; the factor of $1/3$ occurs because, at equilibrium, the conformation is spherical]. If Eq. (2b) is multiplied by $λ/Z$ and the upper convected operator is applied, we obtain the conservative form of the FENE-P dumbbell model¹⁹ in the stress formulation,

Although Eq. (4a) is already much simplified in relation to the original stress formulation of the FENE-P dumbbell model [Eq. (1a)], it is possible to have an even more compact formulation if it is realized that the term involving $aI$ in Eq. (2b) is constant, so it does not contribute to the stress divergence term in the momentum balance. In this case, we re-define stress as $τ≡(ηp/λ)ZA$⁠, and after applying the convected derivative to $(λτ/Z)=ηpA$⁠, it results in

Regarding the conformation formulation of Eq. (2), which corresponds exactly to the original proposal of Bird et al.,^2,3 several variants have been employed in the literature especially that related to turbulence modeling of dilute polymer solutions and the phenomenon of turbulent drag reduction.^31,32 The same equations with the $a$-term substituted by unity with the consequence that neither $A∇$ or $τ$ go to zero at equilibrium, as they should. It may be argued that any such inconsistency is small when $L2$ is large (and $a→1$⁠), but in many applications, lower values of $L2$ are employed (e.g., Purnode and Crochet³³ used $L2=$ 4; also, simulation of concentrated polymer solutions or melts $L2$ is set to small values, see Bird and DeAguiar¹⁶), keeping the original formulation is hardly more involved than setting $a=1$⁠. A possible explanation for this inconsistent simplification is that many authors write the model with a differently normalized conformation tensor, which we call here $B=H⟨QQ⟩/kT$ (since it is the same normalization used for the maximum extensibility parameter $b$⁠). Note that $B$ is not unity at equilibrium [$Beq=(b/(b+5))I$]. The equations become

for the original version without Peterlin's correction and

for the improved version of Bird et al.;² both have $Z=11−TrB/b$⁠.

Clearly, the evolution equation of Eq. (6a) does not have $a$ multiplying $I$⁠, but when the stress is calculated, $a$ is there, and a similar inconsistency arises again. Equations (6b) look more complicated than those for $A$ [Eq. (2)] and are the same as in Bird et al.,² p. 90, except that they denoted $B$ by $α$⁠, the structure tensor.

Other authors (e.g., Vaithianathan and Collins³¹) have normalized the spring function as $Z/Zeq=(L2−3)/(L2−TrA)$⁠, but in this case, $λ$ (and a Weissenberg number based on characteristic velocity U and length scales L_c, i.e., $Wi=λU/Lc$⁠) would need to be adjusted and does not keep the standard molecular definition for the relaxation time, $λ=ς/(4H)$ (where $ς$ is the Stokes drag coefficient on a sphere and H is here the Hookean spring coefficient). However although different from the original FENE-P, such re-scaled FENE-P dumbbell models with $λ′=λ/Zeq$ and $Z′=Z/Zeq$ offer two advantages: (1) the relaxation time becomes $λ′=ψ1,0/2ηp$ and can be determined from the experimentally measured value of the primary normal stress coefficient $ψ1,0$⁠, and (2) the analogy with the sPTT (see Sec. III A) requires just a single condition $ε=1/L2$⁠.

Another meaningful variant of the FENE-P dumbbell model is to substitute the Warner spring force by a “more accurate” expression in the sense of following closer the inverse Langevin function implied by kinetic theory of randomly jointed chains proposed by Cohen,³⁴ namely, $Z=(L2−13TrA)/(L2−TrA)$⁠. In a recent review, Larson and Desai³⁵ are favorable to such an improvement in a FENE-P version for dilute solutions, but if $λ$ is to retain the original definition, $a≡Zeq=(L2−1)/(L2−3)$ should be included in the term multiplied by the identity tensor. Otherwise, the set of equations to be solved is exactly the same as Eqs. (4). Alternatively, it can be employed in the re-scaled FENE-P fashion as explained in the previous paragraph.

Rooted in the network theory for rubber elasticity,^6,7 the PTT model was, in fact, proposed first by Phan-Thien and Tanner earlier than the FENE-P model⁵ [with a function $f(τ)$ linear on the trace of the stress, $τ≡Tr(τ)=τkk$] and afterwards by Phan-Thien³⁶ [$f(τ)$ exponentially increasing with $τ$]. Other variants have followed.³⁷ Following the explanation given by Tanner³⁸ where PTT-like models are contrasted with tube-based models for branched polymers (such as the Pom-Pom of McLeish and Larson³⁹ and the XPP of Verbeeten et al.⁴⁰), the evolution of the conformation tensor in network theories follows:

where $F(A)$ and $H(A)$ are two non-dimensional functions of $A≡Tr(A)$ representing the rates of destruction and creation of network junctions. By assuming that the stress tensor is then given by the simple Kramers'⁴¹ relation (here, we consider only a single mode representation),

then, $λτ∇=ηp(A∇−I∇)=ηp(A∇+γ̇)$⁠, and using Eqs. (7) and (8), Tanner³⁸ obtained

If there is no creation or destruction of junctions as in the case of a dilute polymer solution, then $F=H=1$ and we recover the Oldroyd-B equation (adding a solvent contribution to $τ$⁠). If the rate of destruction is equal to the rate of creation, then $f(τ)≡F=H=1+ε(λ/ηp)Tr(τ)$ and we get the affine (“simplified”) PTT model with the linear function,

Written in this manner, it becomes clear that there is some similarity between this PTT model and the FENE-P dumbbell model of Eq. (4): the effective relaxation time $λ/f$ is now outside the convected derivative (denoting the Z function as f), and the advection of the effective viscosity term $aηp/f$ [the last term in Eq. (4)] is now absent. Finally, if there is no creation of network junctions, $H=1$⁠, and $F(τ)$ takes the expression valid for the XPP of Verbeeten et al.,⁴⁰ then the version of the XPP without the Giesekus⁴² term (⁠$α=0$⁠) is recovered.

The PTT model followed the developments by Tanner²⁵ using the same “delta-function” closure ideas (i.e., the paper where the equations for the FENE-P dumbbell model originally appeared and why Bird proposed calling the FENE-P “Tanner's equation”). Thus, in some sense given by this fundamental link, the PTT and FENE-P dumbbell model may not be expected to be accurate when this closure assumption is invalid. For example, for dumbbell models, although assuming the distribution of dumbbell configurations is highly localized and would appear to be a good assumption for steady-state extension, Keunings³⁷ showed via the direct comparison of stochastic simulations of the FENE model with FENE-P dumbbell model predictions that in transient extension, they exhibit significantly different dynamics of molecular extension and stress. The configuration distribution function for the FENE-P dumbbell model is always Gaussian and, thus, never localized, whatever the flow kinematics.

As a final comment, inspection of Eq. (7) in comparison with Eq. (2), we see that the FENE-P may be seen as a network model with no creation of network junctions (⁠$H=1$⁠). Therefore, we should expect that, in general, the FENE-P dumbbell model and, by the same reason, many tube and reptation-based models, may exhibit much stronger recoil effects as compared with the PTT model. This is, in fact, reported by Tanner³⁸ and will also be observed in the transient flows that we consider here.

Given all the various definitions of the FENE-P dumbbell model in the literature as discussed in Sec. II, in this short section, we explicitly state the form of the models that we are comparing in the remainder of this paper. In particular, we use the simplified linear form of the PTT model (often called the “sPTT” model—a terminology that we will use hereafter—or the affine linear PTT), which can be expressed from Eq. (10) as⁴³ 

where $fPTT=1+λεηpτkk$ is the linear function of the PTT model in which $ε$ is the extensibility parameter. It is the “simplified” Phan-Thien and Tanner model in the sense that only the upper convected derivative is retained (assuming affine deformation), whereas the original model uses the full Gordon Schowalter (GS) derivative⁴⁴ with the parameter ξ, which allows for the slip between the molecular network and the continuum medium. In the sPTT model, ξ is set to be identically zero and the GS derivative simply becomes the upper convected derivative. Note that it is known that the GS derivative can lead to incorrection predictions in some flows,^45,46 and therefore, the slip parameter is often restricted to certain values to avoid such unphysical results as in the model of Johnson and Segalman,⁴⁷ for example.

In contrast, the FENE-P dumbbell model can be expressed as [from Eq. (4) using now f instead of Z for consistency with the sPTT model]

where the symbols have the same meaning for the sPTT model as above, but now, the FENE-P function is $fFP=1+(3a+λτkk/ηp)/L2$ and $a=L2/(L2−3)$ and $L2$ is now called the extensibility parameter. $D/Dt$ still represents the substantial or material derivative [$DfFP/Dt=∂fFP/∂t+(u.∇)fFP$]. Note the the underlined term in Eq. (12) can also be expressed equivalently as $−(λτ+aηpI)fFPD(1/fFP)Dt$⁠.

For any homogeneous steady flow at a constant deformation rate, such as steady-state simple shear flow or planar/uniaxial extension, $DfFP/Dt=0$⁠, and thus, the underlined term on the right-hand side of Eq. (12) will vanish with the FENE-P dumbbell model obtaining a similar expression to Eq. (11) of the sPTT model. For values of $L2≫3$⁠, a tends to one and the function $fFP$⁠, therefore, simplifies to

where now the FENE-P dumbbell and sPTT models can be seen to be identical under the transformation $ε=1/L2$⁠. Thus, in the limit that $L2≫3$⁠, the FENE-P and sPTT material properties (using $ε=1/L2$⁠) should be identical. In fact, as was originally shown by Cruz et al.,¹⁸ there is complete equivalence between the sPTT and FENE-P dumbbell models for this class of flows under the transformation $ε=1/L2, fPTT=fFPa,λPTT=λFPa$ as is immediately apparent when comparing Eq. (4) with Eq. (10) (including the functions). Of course, as L² must remain larger than 3 in the FENE-P dumbbell model, yet no such restriction is required for the ε parameter in the sPTT model, this equivalence only holds for ε < 1/3.

It is useful at this point to look at the response of the various functions, in general conditions (i.e., not simply restricted to homogeneous flows), to try to understand the puzzling fact of how a molecular model whose main assumption is the limitation of its extended length, gives similar results to another model without such a fundamental limitation. Figure 1(a) shows the variation of the “spring” force coefficients $f(x)$ for the FENE-P dumbbell model, including the Cohen version and the theoretical inverse-Langevin function, all being limited by $x=(Q/Q0)2=TrA/L2<1$⁠, and the two PTT versions with linear and exponential functions for which $x=ε(TrA−3)$ may be larger than unity. An explanation for such an unexpected behavior is provided by the stress–conformation relations for the two models: while the limited conformation values in the FENE-P model [Eq. (2b) are multiplied by the force coefficient $f(x)$⁠, which attains large values on approaching the divergence point (⁠$x→1$⁠), in the PTT model [Eq. (8)], the stress is equal to the conformation multiplied by a constant factor for consistency of units. So, when the spring coefficients $f(x)$ are transformed to become functions of the stress $f(τ′)$⁠, where $τ′$ is the trace of the non-dimensional stress normalized by $L2$ or $1/ε$⁠, the results are almost identical as shown in Fig. 1(b). This is the key explanation for the coincidence of results between the FENE-P dumbbell model and the sPTT model in flows of increasing complexity, to be discussed in the next sections.

In this section, we briefly review the equations to be solved and the associated numerical methods used to study our “complex” flows. This encompasses both “start-up” shear, planar extensional, and channel flow and also steady mixed kinematic flows—i.e., encompassing regions of both shear and extension-dominated flows—in five different geometries, including (a) the developing flow in a square duct, (b) the flow in a 4:1 planar contraction, (c) the flow through a cylinder in a channel with 50% blockage, (d) a two-dimensional “cross-slot” geometry, and (e) a fully three-dimensional “Extensional Viscometer-Rheometer-On-a-Chip” (EVROC).⁴⁸ Schematic representations of these geometries are provided in Fig. 2.

The equations solved are conservation of mass, assuming incompressibility, and momentum,

where Re is the Reynolds number, which, in general, is set to be identically zero (the exception is the developing square duct flow, where it is set to be equal to 1). For most cases, the components of the stress tensor are calculated using the so-called log-conformation approach for the sPTT [Eq. (11)] and FENE-P [Eq. (12)] constitutive equations.^49–51 

Rheological constitutive equations in the kernel-conformation approach are generally written in terms of the conformation tensor introduced in Sec. II, $A$⁠, which is a variance–covariance, symmetric positive definite tensor (SPD)⁵¹ and in the non-dimensional form [from Eq. (8)], we have

for the sPTT model and

for the FENE-P model [from Eq. (2b)], where the functions from the FENE-P (⁠$fFP$⁠) and sPTT (⁠$fPTT$⁠) models in terms of the conformation tensor are

To solve Eqs. (14)–(18), three different codes have been used (depending on the problem studied). For the developing flow in a square duct and the EVROC geometry, we use the in house finite-volume code as described by Zografos et al.⁵² and full details are provided there. For the start-up channel flow problem, the flow is generated by a sudden application of a pressure gradient, and the equations to be solved are given in the supplementary material of Alves et al.⁵³ or Oliveira⁵⁴ [we highlight that for the FENE-P dumbbell model, one needs to solve the polymer stress equations for the lateral normal stresses (⁠$τyy=τzz$⁠) despite the two-dimensional nature of the flow field]. For the remaining geometries, and the homogeneous “start-up” simulations, we use the open source RheoTool package,⁴⁹ where a Gauss scheme is used to discretize the gradient, divergence, and Laplacian terms. In these cases, a zero gradient boundary condition for the stress and velocity components is applied at the outlets to simulate fully developed conditions. At the walls, values of the stress components are calculated using an extrapolation method as suggested by Pimenta and Alves,⁴⁹ and no-slip is assumed for the velocities.

For the 4:1 contraction, the flow through the cylinder and the cross slot the geometries used in this work are identical to the ones used in the tutorial section of the RheoTool package,⁴⁹ and for the purpose of brevity, information regarding the mesh and boundary condition is not unnecessarily repeated here. Mesh details for the square duct (30 × 30 cells in the cross plane) are given in Table I. The boundary conditions and EVROC geometry mesh (M1) are the same as in Zografos et al.⁵² 

As discussed earlier the FENE-P dumbbell and sPTT models can be seen to be identical under the transformation $ε=1/L2$ for simple shear and extensional flows. Thus, in the limit that $L2≫3$⁠, the FENE-P and sPTT material properties (using $ε=1/L2$⁠) should be identical. [For larger(smaller) values of $ε(L2)$⁠, complete equivalence¹⁸ between the sPTT and FENE-P models occurs when $ε=1/L2,fPTT=fFP/a,λPTT=λFP/a$⁠.] We show this in Fig. 3 for a range of conditions, including solvent viscosity contributions (although as this is a linear addition, this precise value of solvent contribution is insignificant in demonstrating this equivalence). Exact analytical expressions for the sPTT material functions are provided in Shogin²² and are not unnecessarily repeated here.

Once a model's shear viscosity response is known, this information is then sufficient⁵⁶ to obtain analytical (or semi-analytical) solutions to any steady fully developed simple-shearing flow such as Couette flow or laminar flow through axisymmetric pipes, two-dimensional channels, or concentric annuli. Thus, solutions obtained for the sPTT model in such flows, but not yet independently for the FENE-P model, do not need to be separately derived. Examples for the sPTT model include annular flow solutions due to cylinder rotation in a concentric annulus⁵⁷ or channel flows with slip,^58,59 electro-osmotic Poiseuille flows,^60–62 or thin-film flows.⁶³ For the FENE-P model, there are boundary-layer and free-shear layer approximations,^64–66 which should be equally valid for the sPTT model.

In Fig. 4, we also shown that this equivalence holds in the planar extensional flow under the same conditions.

For the start-up flows considered in this section, the strain rates (i.e., $γ̇$ and $ε̇$⁠, respectively, for shear and extensional flows) can be written as

where $H(t)$ stands for the Heaviside step-function with $H(t)=1$ for $t>0$ and $H(t)=0$ at $t=0$⁠. Such flows have been studied previously using both numerical^19,67,68 and analytical approaches^21,22 although direct comparisons between the sPTT and FENE-P dumbbell models have not been previously made to the best of our knowledge.

The results presented in Figs. 5 and 6 show the planar extensional and shear viscosity response of the sPTT and FENE-P dumbbell models under start-up flow conditions that are numerically obtained using the rheoTestFoam solver introduced in the Rheotool package.⁴⁹ (We confirmed that for the selected parameter values, the results are in good agreement with the analytical solution of Ref. 22 for the sPTT and numerical results of Ref. 69 for the FENE-P dumbbell model). The results show that as the Weissenberg number is increased, the long-time (i.e., steady-state) extensional and shear viscosities exhibit larger and smaller values than the zero-shear rate viscosity, respectively, which is related to the tension-thickening and shear-thinning behavior of these models as already shown in Figs. 3 and 4. As already discussed for the steady-state results, the long-time (steady-state) values of both shear and extensional viscosity are identical in both the sPTT and FENE-P dumbbell models under the transformation $ε=1/L2$(in the limit of $L2≫3$⁠). As can be seen in Figs. 5 and 6, the equivalence of the two models which was observed under steady-state conditions does not generally hold in transient flow conditions. In particular, as you move outside of the dilute regime condition [i.e., $(1−β)$ must remain small] and especially for larger Wi values, there are marked differences in the transient response of the two models. In the limit of $L2≫3$⁠, due to the presence of the extra material derivative in the FENE-P dumbbell model [the underlined terms in Eq. (12)], the deviation of the two models from the Oldroyd-B model is not identical under transient flow conditions. For the planar extensional flow, it can be seen that around Wi ∼ 1 differences in the transient response between the two models can be seen regardless of the value of $β$⁠. In contrast for the shear flow, which is known to be a weaker form of deformation than extension,⁷⁰ the transient response is seen to be approximately the same between the two models until Wi ∼ 2.

Our results suggests that for extension during the start-up process, increasing the Weissenberg number, the time required to reach steady-state conditions is reduced (notice that time is normalized by the extension or shear rate) in both models, while the FENE-P dumbbell model exhibits very sharp transitions at the highest Wi studied. For $Wi>2$⁠, the extensional viscosity exhibits an almost step-like change at a critical time $ε̇0 t≈3$⁠. This behavior looks more pronounced as the viscosity ratio becomes smaller because of the logarithmic scaling used in the plot but is essentially the same for the two solvent viscosity ratio cases shown. This rather unphysical behavior of the FENE-P dumbbell model in transient extension has been observed previously and modifications suggested (e.g., FENE-PB model⁶⁹ or FENE-M,¹⁹ which is equal to the sPTT model when the two non-linear spring force functions coincide).

Further analysis of this behavior is possible on the basis of the governing constitutive equations for planar extension for the FENE-P dumbbell model, and the one of importance here is that giving the evolution of the axial normal component of the conformation tensor $Axx$ (since $Ayy≪Axx$⁠). Full details of this analysis are given in the  Appendix. From Fig. 5, it is clear that for sufficiently high $Wi≡λε̇$⁠, where $ε̇=∂u/∂x$ is the elongational rate (in practice “high” is when $Wi>2$⁠), stretching becomes exponential, that is the equation simplifies to

and by assuming that the $f(Akk)≅f(Axx)$ function remains approximately constant at its minimum value, $feq=a$ gives

and then $Axx(t)$ evolves very quickly up to the limiting value imposed by the spring stiffness function [this being $Axx,max=(L2(2Wi−1)−a)/(2Wi−13)$ with the 1/3 becoming 0 if the Cohen spring is not considered; $ηE,max=2Axx,max$]. Note that the timescale for normal stress and consequently elongational viscosity growth is $≈1/(2λε̇−1)=0.33$ for $Wi=2$ and $≈0.05$ for $Wi=10$ (in units of $λ$⁠). For lower Weissenberg numbers $Wi≤2$⁠, the complete equation for $Axx$ needs to be considered, which is approximated to first order as

The solution is now an hyperbolic tangent function (⁠$m=2Wi−1$⁠; $m′=m2+4a/L2≈m$⁠; $t1$ is the time for $Axx/L2$ to be half $m$ and is defined in the  Appendix),

which is very closely followed by the numerical simulations up to the point where $Axx$ attains its maximum value (at $t=ts$⁠). The analysis also allowed us to estimate this time for maximum molecular stretch; for large $Wi$⁠, it is $ε=ε̇ts=log L$⁠. Such a behavior is clearly seen in the curves of the related quantity of extensional viscosity vs $ε̇t$ for $Wi≤2$ in Fig. 5 [in the late time steady-state, the planar extensional viscosity is simply $2Axx≈2L2(1−1/2Wi)$]. As $Wi$ gets larger, the maximum $ηe$ is reached when the Hencky strain attains $ε=ε̇ts=3.1$ for $L2=500$ of Fig. 5.

In contrast to extensional flows, the start-up shearing-flow reaches steady-state conditions faster at lower Weissenberg numbers. In marked contrast to the extensional data, which remained monotonic, in Fig. 6, there are noticeable “overshoots” of the transient viscosity at higher Weissenberg numbers (these also show up in the transient normal stress data, which is not shown for conciseness). Such phenomena have been studied within the framework of Oldroyd's 8-constant framework⁵⁸ although neither the sPTT model or the FENE-P model falls within this hierarchy of models as they contain terms involving τ_kk. These overshoots (and undershoots) at high values of the Weissenberg numbers occur in both models but these oscillations are shown to be significantly more pronounced in the FENE-P dumbbell model. As previously noted,¹⁹ the very large overshoots in shear of the FENE-P model—where the overshoots exceed the linear viscoelastic envelope³—are not observed in experiments for polymeric fluids.⁷¹ 

This problem represents an extension of the stress growth upon inception of the simple shear flow of Sec. V B to the non-homogeneous velocity-gradient case where the flow is generated by a sudden application of a pressure gradient. It requires solving the equation of motion and has the advantage, compared to the complex flows of Sec. V D, of isolating time-dependent elastic effects. A constant pressure gradient $P≡−dp/dx$ is applied at $t=0$ to the fluid initially at rest. In non-dimensional terms, we choose $P=3$ so that a constant-viscosity model will have an average cross-section velocity of unity (⁠$U=1$⁠) when a fully developed situation is attained at large times. This means that in the results to be shown, lengths are scaled with the channel half-width $H$⁠, velocities with the characteristic velocity $Uc=PH2/3η0$⁠, and time with a diffusion timescale $tD=(ρH2/η0)$⁠. Non-dimensional parameters defining the solution are the elasticity number, $E=λ/tD=λη0/ρH2$⁠, which we set at the moderate value of $E=5$⁠, the Weissenberg number $Wi=λUm/H$ (⁠$Um$ is the average velocity at steady state), which is indirectly set by the choice of $P$⁠, and the solvent viscosity ratio, $β=ηs/η0$⁠.

Figure 7 shows the evolution of the centerline velocity $u0=u(y=0,t)$ for the FENE-P dumbbell and the sPTT models using two typical values of the extensibility parameter $L2=100$ (thus, $ε=0.01$ for sPTT) and $L2=5000$ (⁠$ε=0.0002$⁠). Due to shear thinning, for the same applied pressure gradient, the flow rate (and average velocity, $Um=Q/2H$⁠) will be larger for the two viscoelastic models, especially as $L2$ gets smaller (and $ε$ larger). Hence, we see that by $t=40$⁠, the flow has attained a fully developed condition with $u0=1.5$ for the Oldroyd-B fluid, while the two viscoelastic models have exactly the same fully developed solution with $u0=1.559$ for $L2=5000$⁠, still showing a small shear thinning effect (normalized with the average velocity $u0/Um=1.489$⁠), and $u0=3.125$ for $L2=100$ (⁠$u0/Um=1.384$⁠). During the transient regime, the results for the FENE-P and sPTT models differ considerably, especially for low $L2$ and high $ε$⁠, and the differences last for about 3–4 relaxation times (up to t = 15 to 20). We may also do calculations with differing pressure gradients so that the final average velocity at steady state will be the same (⁠$Um=1$ in dimensionless form). From the analytical solution, we find that $P=1.836 02$ for $L2=100$ and $P=2.875 75$ for $L2=5000$ (keeping $E=5$⁠). Note that P should rise up to $P=3$ as the effect of shear thinning becomes negligible. Figure 8 shows the results of these calculations, which correspond to $E=Wi=5$⁠, and various imposed pressure gradients.

The evolution of the centerline velocity is rather similar to that seen previously (Fig. 7), even more so if Fig. 7 is represented as $u0/Um$⁠. For $L2=5000$⁠, the FENE-P and sPTT results almost coincide with each other during the whole transient evolution. For $L2=100$⁠, there is coincidence only after about $t>≈15$ (after three relaxation times), and the velocity undershoot seen for the FENE-P dumbbell model after the initial almost inertial velocity increase (and stretch of the fluid molecules) is much more pronounced than for the sPTT model. If we identify that velocity decreases as a manifestation of recoil, Fig. 8 shows that recoil for the FENE-P model is more intense than recoil for the sPTT model.

As most complex flows involving extensional flow kinematics only involve transient stretching due to finite residence time/Hencky strain (with the noticeable exception of cross-slot type flows containing internal stagnation points that couple finite strain rate and zero velocity⁷²), the results in this section and the Sec. V B would suggest that non-linear simulations of the sPTT and FENE-P dumbbell models in complex geometries may exhibit significant differences. We now test this hypothesis in a range of such geometries in Sec. V D.

In Subsections V A–V C, we have shown that, although in steady-state homogeneous flows both models predict identical behavior at small extensibility parameters under the transformation $ε=1/L2$⁠, beyond a critical Weissenberg number, they exhibit different stress growths under otherwise identical transient “start-up” tests in both homogeneous shear and extension and non-homogeneous start-up of the channel flow. Start-up simulations for both shear and extensional regimes show that agreement between the models only holds under transient conditions until Wi ∼ 2 for the conditions highlighted.

In this section, we extend our analysis to more geometrically complicated flows by undertaking a series of comprehensive simulations in complex geometries. We restrict our analysis here to flows that remain Eulerian steady (so typically, Wi based on characteristic length and velocity scales remains on the order of 1). Of course, these flows are Lagrangian unsteady and that is why we refer to them as “complex.” The flows selected are those which have previously been proposed as “benchmark” geometries as they are known to be challenging and contain significant regimes of both shear and extension-dominated kinematics.⁵⁵ 

In this section, the behavior of the FENE-P dumbbell and the sPTT models is investigated for a flow within a duct of square cross section. For all cases, a uniform velocity profile is applied at the inlet of the geometry (magnitude U), and thus, the duct is made long enough so that the flow is able to reach to a fully developed state. This length is L = 60H, where H is the width (and the depth) of the duct. The Weissenberg number for this configuration is then defined as Wi = λU/H. Although this flow is seemingly complex, this is only true in the developing region of the flow as, once it is fully developed, this flow is basically “viscometric”⁴⁶ and the response is completely determined via the shear viscosity (the normal stress $τxx$ playing no role in the momentum equations). This holds because neither the sPTT model nor the FENE-P dumbbell model exhibits a non-zero second normal-stress difference in steady simple shearing.⁷³ Thus, developing steady flow in a square duct represents an interesting analog of the transient start-up flows discussed in Sec. V D in that the flow develops spatially and is a mixture of both shear and extension (as the flow is uniform at th inlet, it must be decelerating simple shearing near the walls and acceleration/purely extensional along the centerline). Finally, just as the start-up flows at long times must approach the steady-state values, in this case, the flow develops spatially until it becomes fully developed (i.e., there is no further spatial development), and then, the flow is simple shear and also equal to its steady-state viscosity (at the local shear rate magnitude). However, in contrast to the homogeneous flow fields discussed in Sec. V A, the shear rate (and hence viscosity) is non-homogeneous and increases from zero at the centerlines to larger values at the walls.

In Fig. 9, the normalized streamwise velocity profiles along the centerline [Fig. 9(a)] and along the transverse direction on the centreplane [Fig. 9(b); z = 0] for the cases of ε = 0.01 and L² = 100, β = 1/9 are shown. It can be seen that for all the investigated Weissenberg numbers, when the flow becomes fully developed, both models predict very similar maximum velocities and identical transverse profiles (maximum velocities agreeing to within better than 0.16% in all cases). More interestingly, the developing flow near the inlet—as evidenced by the variation of the centerline velocity in Fig. 9(a)—shows that there are marked differences between the models. Although both model data sets exhibit significant overshoots close to the inlet (reaching ∼ 2.1U), the sPTT data decay monotonically back to the fully developed values, which themselves decrease with the increasing Weissenberg number (as a consequence of shear thinning). In marked contrast, the FENE-P data at higher Weissenberg number first overshoots then undershoots the fully developed values. As the flow along the centerline is pure-extension (the flow is accelerating, and due to symmetry, there must be no shear), we see over-and under-shoots in velocity, which were absent in the equivalent start-up extensional flow. These interesting features must be a consequence of the non-homogeneous nature of the extensional deformation along the square-duct centerline. The behavior within the same configuration when L²= 5000 for the FENE-P model and ε = 0.0002 for the sPTT model, and the obtained profiles, as shown in Figs. 9(c) and 9(d), are basically identical and Oldroyd-B like.

In this subsection, we present an investigation of the viscoelastic flow in a planar sudden contraction of contraction ratio 4 (i.e., the downstream channel height is simply a quarter of the upstream channel height). This problem is a very well-studied “benchmark” in non-Newtonian fluid mechanics,⁴³ and therefore, we fix the solvent-to-total viscosity ratio at β=1/9 in order to be consistent with those conditions. We test both the FENE-P and sPTT constitutive models to compare the difference between the results in the limit that $L2$ is large and ε is small. Although both models have been directly compared in this geometry previously,^74,75 these studies have not focused on this same parameter range.

The 4:1 contraction is a flow in which the upstream flow, once fully developed away from the inlet, is simple shearing and then the acceleration bought by the presence of the contraction gives rise to locally more complex flow, mixing regions of shear, rotation and, along the centerline, pure (planar) extension (see, e.g., contour plots of the flow-type parameter in Mompean et al.⁷⁶). We note that the nature of the flow along the centerline is always transient extensional flow as fluid elements are being advected at the same time as being stretched, such that they do not accumulate a significant amount of Hencky strain. Far downstream of the contraction, the flow once again becomes simple shearing. From a fundamental point of view, the viscoelastic fluid flow through such ducts with the abrupt change of cross section, either contractions (as studied here), expansions,⁷⁷ or contraction/expansions,⁷⁸ is important as they highlight many unusual phenomena brought about by elasticity. These phenomena include complex recirculation patterns, vortex enhancement or suppression, the possibility of the unsteady flow due to elastic instabilities, complex stress behavior near geometrical singular points, and enhanced pressure drops significantly above the equivalent creeping flow Newtonian result.⁷⁹ 

In Fig. 10, we show streamlines superimposed on top of contours of the trace of the polymeric stress tensor (to give an unambiguous measure of the “elasticity”) which show that, qualitatively, the flow at Wi = 3 is very similar between the two models. To highlight the small differences, in Fig. 10(c), we show a contour plot of the differences between these trace fields, which do highlight some deviations that are concentrated around the sharp corners. More quantitatively, the variation of the corner and lip vortex sizes represented as $Xr$ and $Xl$⁠, respectively, are shown in Fig. 11 (these quantities are defined in Alves et al.⁴³). By increasing the Weissenberg number (defined here based on the downstream channel height and average velocity here), the corner vortices reduce, whereas the lip vortices increase in size. In general, the difference between the sPTT and FENE-P dumbbell models is small for the selected conditions especially for the length of the corner vortices. On the other hand, the difference between the sizes predicted by sPTT and FENE-P dumbbell models for the lip vortex seems to be slightly smaller when $ε=0.01$⁠, which could be related to the singularity appearing near the re-entrant corner.⁴³ 

The flow through a cylinder confined in a planar channel is another classical benchmark problem in non-Newtonian fluid mechanics.^53,80,81 Although there is no geometric singularity in the flow, unlike the sharp corners present in the planar contraction flow, for example, the interior stagnation points, especially that at the rear of the cylinder, may drive the growth of very large polymeric stresses downstream of the cylinder especially along the centerline.^82–84 Thus, this problem is particularly challenging from a numerical simulation viewpoint and should be a case where differences may be expected to occur between the two models due to this strongly Lagrangian unsteady flow along this line. Representative streamlines superimposed on contours of the trace of the polymeric stress shown in Fig. 12, once again highlight that qualitatively the two models—here shown at Wi = 0.5 (note here, for consistency with the literature, Wi is based on the radius of the cylinder and upstream velocity)—are in very close agreement (also highlighting that, under these conditions, the largest stresses appear in the regions of shear above and below the cylinder in the gap where the flow must accelerate). Figure 12(c) shows that differences between the two fields are also concentrated more toward the poles of the cylinder rather than in the birefringent strand downstream of the rearward stagnation point.

In Fig. 13, the variation of the drag force with the Weissenberg number for three different viscosity ratio parameters of $β=0, 1/9$⁠, and 0.59 are reported for two different extensibility parameters of $ε=1/L2=$ 0.01 and 0.002. As can be seen, the difference between the predicted drag coefficient for both models are small and the agreement between the two models gets slightly better as we approach the dilute regime limit (⁠$β→1$⁠). Following this strong agreement between the two models for the prediction of the drag coefficient, we show the variation of the streamwise normal stress $τxx$ along the symmetry line downstream of the cylinder in Fig. 14. The variation of the normal stress with $ε=1/L2=$ 0.002 is almost negligible between the two models. At a higher value of $ε=1/L2=$ 0.01, the difference between the two models becomes slightly larger but again as you move toward the more dilute regimes, the difference becomes smaller.

The cross-slot geometry is a stagnation point flow, which allows large strains to develop in a well-controlled flow field. It can, therefore, be used for extensional rheometry measurements⁸⁵ or once instability arises, as a mixing device.⁸⁶ Indeed, it has been proposed as a benchmark case precisely because the base steady symmetric flow may undergo a bifurcation to a steady asymmetric configuration at well-defined conditions.⁸⁷ Despite the presence of the extensional flow at the stagnation point and the concomitant large extensional stresses, which may accrue at this location, it has recently been shown that the steady symmetry-breaking instability is, most likely, related to the (transient) shear-dominated flow around the corners.^72,88 The critical value of the Weissenberg number where this steady asymmetry occurs, thus providing a subtle differentiator of constitutive equations (see, e.g., the FENE models compared in Ref. 89). In Fig. 15, we show the variation of the asymmetry parameter (A_p) (see the definition provided in Ref. 90) with the Weissenberg number for two different values of $ε=1/L2=$ 0.01 and 0.002. Although at the larger value of $ε$ (smaller $L2$⁠), there are some noticeable differences being on the order of 10% under the critical conditions; as $ε$ is reduced, the two models essentially bifurcate at the same critical value of Wi (∼0.37).

Microfluidic contraction devices have been proposed for extensional rheometry measurements, in particular, as a useful method for determining the extensional viscosity of low elasticity solutions. The first commercially available “extensional viscometer-rheometer-on-a-chip” (e-VROC™)⁴⁸ developed by Rheosense is a hyperbolically shaped contraction/expansion geometry, which incorporates pressure-drop measurement capabilities. We have previously studied this geometry numerically,⁵² including comparisons between the sPTT and FENE-P dumbbell models (see, e.g., Fig. 13 in Zografos et al.⁵²).

Figure 16 shows typical results through the EVROC geometry. At small L² (=100), the FENE-P dumbbell model is shown to produce a very large velocity overshoot and concomitantly increase the strain rate [Figs. 16(a) and 16(b)]. The pressure drop for these model conditions is less than or approximately equal to the equivalent for a Newtonian fluid [Fig. 16(c)], and it appears that shear thinning dominates the pressure drop over any extensional or transient effects, especially when β = 1/9. Small differences in the pressure drop between the two models are apparent at the highest values of Wi that could be achieved in these cases (∼4). At the same value of solvent to total viscosity (β=1/9), but at significantly higher L²(= 5000) (velocity and strain rate profiles not shown for conciseness), the two models behave very similar in terms of both velocity profiles and pressure [Fig. 16(c)]. At higher levels of solvent viscosity (β = 0.95), the pressure drop for the sPTT model essentially remains the same as the Newtonian result for all Weissenberg numbers. In contrast, the FENE-P dumbbell model, beyond a Wi ∼ 10, starts to show pressure-drop enhancements. As the flow along the centerline is basically the same between the two models, the differences observed in the dissipation here must be due to the shear flow outside of this region.

Others have previously discussed similarities between the sPTT and FENE-P dumbbell models in homogeneous flows,¹⁸ between the PTT model and the extended pom pom (XPP) model,⁹¹ and also between the modified FENE-P equation and the sPTT model¹⁹ but a comparison across the range of flows that we have discussed in the current paper has not previously been undertaken.

Although the response of the two models in steady viscometric flows is essentially identical for $ε=1/L2$ when $L2≫3$⁠, this does not necessarily mean results in more complex flows—where the underlined term in Eq. (12) is not identically zero—will necessarily be the same. To investigate how much the models may differ in more general flow conditions, we investigated the response of the two models when the flow is “complex” in a number of different definitions: first, when the applied deformation field is homogeneous in space but transient in time (so-called “start-up” shear and planar extensional flow), as an intermediate case, the start-up of the channel flow and then for a range of “complex” flows, which, although Eulerian steady, are unsteady in a Lagrangian sense. In such flows, the effects of unsteadiness or non-homogeneous flow/Lagrangian unsteadiness may result in differences between the models and indeed cause the functions $fPTT$ and $fFP$ to differ, for example, if complexity in the flow/flow history causes $τkk$ to differ between the two models. Despite these differences in the two models, we find fairly good agreement between the two models, especially in the complex flows where we restricted our simulations to the (Eulerian) steady regime. The equivalence in the spring functions between the two models help to explain these similarities, despite the very different molecular origin of the two models.

R.J.P. would like to thank the EPSRC for the award of a Fellowship under Grant No. EP/M025187/1. P.J.O acknowledges funding from FCT-Fundação para a Ciência e a Tecnologia through project UIDB/00151/2020. Professor Gareth McKinley first suggested looking at differences between the two models in start-up flows (as described in Sec. V B) and this is gratefully acknowledged.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

To better understand the behavior of the FENE-P dumbbell model observed in the numerical “start-up” simulations of transient extensional flow, in this appendix, we analyze the underlying equations based on the conformation tensor [Eq. (2a)]. For transient uniaxial or planar extension, the evolution of the primary normal component of $A$ is governed by

The Warner function [Eq. (3)] is $f(x)=1/(1−x)$⁠, which is a monotonically increasing function of the square of stretch $x=Akk/L2$⁠. If we assume that $Ayy≪Axx$ (verified by simulations provided $Wi≡λε̇≥≈0.1$⁠), then $Akk≅Axx$ or $x≅Axx/L2$⁠, and Eq. (A1) written for $x=Axx/L2$ becomes

When $Wi≫1$ (in practice, as we see in Fig. 5, $Wi>2$⁠), the uncoil of the molecule is exponentially fast and it may be assumed that $f$ remains approximately constant at its equilibrium value $feq=a=L2/(L2−3)$⁠, giving the largest value for the parameter $m$ defined in Eq (A2). With this assumption and neglecting the effect of the last term $a/L2$⁠, the solution, starting from equilibrium (initial condition: $Axx=1$ and $x=x0=1/L2$⁠, for $t=0$⁠), is

Figure 17 shows that the solution provided by Eq. (A3) captures very well the variation observed in the full numerical simulations. At these large values of $Wi$ and $L2$⁠, the elongational viscosity $ηE$ that was shown in Fig. 4 basically follows $Axx$ (⁠$ηE≈aAxx/Wi$⁠) or can be determined more precisely as $ηE≅fAxx/Wi=L2x/[(1−x)Wi]$ with $x$ from Eq. (A3).

The time instant at which the molecule is totally stretched (⁠$ts$⁠) may be evaluated by intersecting the exponential growth of Eq. (A3) with the final limiting value of $Axx$ [obtained by solving the quadratic equation that arises from Eq. (A1) when $∂Axx/∂t=0$⁠, approximately $Axx,∞=(L2(2λε̇−1)−a)/2λε̇$],

For example, for L²=500 and $Wi=λε̇=20$⁠, Eq. (A4) gives $ε̇ts=3.17$⁠. At large $Wi$⁠, we have the following limiting values:

where $ε=∫ε̇dt$ is here the accumulated strain (not to be confused with the PTT parameter). Hence, the case $L2=500$ has a maximum stretch or Hencky strain of $εmax=ε̇ts,∞=3.107$⁠. This explains the behavior observed in Fig. 5. The parameter $L2$ changes the maximum attainable $Axx$ but does not otherwise change the described behavior, as shown in Fig. 17 with $L2=100$⁠. When the Weissenberg number is not large, say $Wi≤2$⁠, the constant term $a′=a/L2$ in Eq. (A2) cannot be discarded and, considering that $x≡Axx/L2$ is initially small, we obtain

whose solution is (now, $m=2Wi−1$⁠)

with

The complete approximated solution to the normal conformation component is finally obtained by limiting the value of Eq. (A7) with the large-time, steady-state solution (⁠$Axx,max$ used above),

Figure 18 illustrates the improvement achieved by Eqs. (A8) and (A9) [denoted 1+tanh(m′t) in the figure] compared to the previous exponential approximation of Eq. (A3) [denoted exp((2Wi-a)t)]. It is also apparent from Fig. 18 that for Wi > 2, Eq. (A3) becomes adequate.

For the sPTT model with the linear function, the solution for $x=εAxx$ is also given by Eq. (A7) but without the need of the limitation imposed by Eq. (A9); the new meaning of the parameters are $m=2Wi−1+83ε$⁠, $a′=ε(1−53ε)$⁠, and $Axx,∞=m/ε$⁠.