Viscoelastic laminar drag bounds in pipe flow

Elasticity effects have long been known to affect transition to turbulence and drag in fluid flows, whether these effects are induced by fluid–structure interaction^1–3 or by viscoelastic rheological aspects associated with non-Newtonian fluids.^4–8 For instance, the addition of low concentrations of long-chain polymers generates an astounding 80% drag reduction (DR) in turbulent regimes,^9,10 which has significant implications for practical applications. At relatively low Reynolds numbers, elasto-inertial effects set in,^11,12 which are also observed in the instabilities of flow bounded by compliant walls.^1–3 

When considering turbulent dilute polymeric solutions, the exhibited DR has been attributed to the interplay between flow turbulence and elasticity of the polymers in the near-wall region.^9,10,13,14 DR is most prominent when the time scale of the polymer elastic dynamics—known to be dependent on the number and length of monomers making the polymer—is of the same order or higher than that of the turbulent fluid flow.¹⁵ In addition, such fluids exhibit a maximum drag asymptote with respect to the polymer concentration,^15–20 which is suggested to be associated with elasto-inertial instability.^11,12 A phenomenon of reverse transition has even been uncovered in such flow recently.^21,22 The energy cascade is also different from that of the Newtonian counterpart²³ since some energy is rerouted to the polymer stretching dynamics, thereby reducing the formation of the smallest eddies and the associated viscous energy dissipation.

However, this DR phenomenology is absent in steady laminar flows with polymers^9,16 due to the absence of small time scale in the flow dynamics. Here, we show analytically that the laminar drag of a FENE-P (finitely extensible nonlinear elastic model under Peterlin approximation) fluid (the solution) in a cylindrical pipe exhibits a set of lower and upper bounds. Specifically, the laminar drag is lower than that of a Newtonian fluid with a viscosity matching the total viscosity of the solution, while being higher than that of the pure solvent. This previously unreported phenomenon is due to the effective viscosity of the solution being bounded by the limits of that of two Newtonian flows: (i) the flow of pure solvent and (ii) that of a viscosity-matched fluid.

Finally, we relate the radial component of the stress tensor to the effective viscosity, which has observational support from direct numerical simulations of turbulent flows.²⁴ Based on this result, we give an interpretation of shear-thinning, which is well-known to be directly correlated with the axial elongation of polymer molecules. We argue that the perceived effect of shear-thinning is due to a sequence of constituent fundamental physical phenomena, such as elasticity and force balancing wrapped together.

The flow dynamics is characterized by three nondimensional parameters: (1) Re = ρU_cR/μ, (2) β = μ_s/μ, and (3) Wi = λU_c/R, with μ (resp. μ_s) being the total viscosity (resp. the solvent viscosity), U_c the centerline velocity in the absence of polymer, R the pipe radius, and λ the elastic relaxation time. A steady flow of dilute polymers modeled as dumbbells is governed by²⁵ 

where τ = (fc − I)/Wi is the elastic stress of the polymers of maximum extensibility L. Under the FENE-P model, the Peterlin function takes the form f = (L² − 3)/(L² − tr(c)). The conformation tensor has, for entries, $cij=⟨R̃iR̃j⟩$⁠, where $R̃i$ is the end-to-end vector of a polymer molecule and I denotes the identity matrix. For a formulation of the FENE-P model alternative to Eqs. (1) and (2) with an extended set of parameters and a different set of unknowns, refer the study of Bird et al.,²⁶ which has been solved by Cruz et al.²⁷ for laminar profiles in the pipe and channel. (For the profiles for inviscid solvent, refer the study of Oliveira²⁸). However, the FENE-P model is widely known to the community studying turbulence and transition as it appears in Eqs. (1) and (2) with the Peterlin function f as mentioned. We therefore solve these equations below and use the resulting analytic solutions to study the bounds for velocity and friction. The laminar profiles derived here will also enable the study of the transition to turbulence by perturbing Eqs. (1) and (2). Hitherto, the absence of a solution to Eqs. (1) and (2) has resulted in that transition studies have only been conducted under the Oldroyd-B model.⁷ 

For steady laminar flows, p = P(x), c = C(r), and u = [U(r), 0, 0] using cylindrical coordinates (x, r, θ), which are in axial (ê₁), radial (ê₂), and azimuthal (ê₃) directions, respectively. The components of $∇$·τ in Eq. (1) can be found in the study of Bird et al.,²⁹ while Eq. (2) becomes

where F(r) = (L² − 3)/(L² − tr(C)) and the subscript r denotes the derivative along the radial direction. Equation (3) gives the non-zero elements of C,

Equation (5) implies that the transverse components of τ follow τ₂₂ = τ₃₃ = 0, which is in agreement with Ref. 27. Substituting Eqs. (4)–(6) in F yields

Since the driving pressure gradient in a pipe is given by

Eq. (1) can be written as $βrUrr+(1−β)rUr/Fr=−4r,$ which under the condition U_r(0) = 0 gives

Substituting this expression for U_r into Eq. (7), we have

Given that Eq. (7) implies F − 1 ≥ 0, a solution of Eq. (10) reads

where $ζ1=a+a2−1$⁠, $ζ2=a−a2−1$⁠, and a = 1 + 108βWi²r²/L² with the positive square-root and real cubic-root implied. Eventually, Eq. (9) gives the velocity profile

where $U(r)$ is the Heaviside step function. The set (4)–(6), along with Eqs. (11) and (12), gives the complete FENE-P steady laminar pipe flow solution that we use to study the laminar drag. Their behaviors close to the pipe center, i.e.,

show that they exhibit even or odd symmetries with respect to r. These parities can assist when analyzing the symmetries of small perturbations close to the wall as performed for the case of Newtonian flows³⁰ and when deploying them in the numerical calculations.^31,32

Figure 1(a) shows the velocity profiles in the Newtonian case (β = 1) and three non-Newtonian cases (β = 0.8 and increasing Wi). They reveal that under the same driving pressure gradient, the volume flow rate $q̇=2π∫01U(r)rdr$ increases when β drops from 1 to 0.8, thereby implying an apparent DR in the laminar regime. This increase in $q̇$ is further amplified when considering increasingly large Wi. However, this should not be interpreted that the addition of polymer reduces drag compared to the pure solvent Newtonian case.

To better understand that, let us recall that the pressure gradient (8) is defined based on the Reynolds number and therefore the total viscosity μ. For a given fixed value of dP/dx, the total viscosity has to be the same for all values of β. Hence, when β is altered, the viscosity of the solvent changes for the profiles shown in Fig. 1(a). This implies that the various profiles for different β values correspond to different solvents and that the non-Newtonian flow effectively exhibits less drag when compared to a Newtonian fluid—other than the solvent—having the same total viscosity as the one at the particular β value considered.

We now show that the drag of a Newtonian flow of the pure solvent is lower than that of the non-Newtonian solution after adding polymers. In the non-Newtonian case, the dimensional pressure gradient is given by dP^*/dx^* = −4μU_c/R², where the superscript “*” indicates the dimensional value of a given variable. If the same dimensional pressure gradient is applied on the pure solvent, the centerline velocity increases by a factor of 1/β since dP^*/dx^* = −4μ_s[U_c/β]/R². Hence, the dimensional velocity profile of the pure solvent is $Us*=Uc(1−r2)/β$⁠, and its non-dimensional counterpart U_s(r) is shown in Fig. 1(b) together with the profile U(r) of the polymeric solution. As can be clearly observed, the pure solvent experiences less drag. To study a range of Weissenberg numbers Wi, we compare the volume flow rate of the polymeric solution, $q̇=2π∫01rU(r)dr$⁠, against that of the pure solvent, $q̇s=π/(2β)$⁠. As shown in Fig. 1(c), there is always a drag enhancement for all values of Wi, and this result holds when L is varied.

In the limit of Wi → ∞ or L → 0, we find F → ∞ from Eq. (7). Hence, considering the same limit in Eq. (9), we obtain U_r → −2r/β, which is the derivative of (1 − r²)/β. These results therefore provide the bounds for the velocity profile as

for all values of Wi and L. Physically, this implies that the polymeric solution experiences less drag than the Newtonian flow of a fluid having a viscosity identical to the total viscosity of the former, but higher drag than that of the pure solvent under the same driving pressure gradient. This can be made more clear from the behaviors of their respective friction factors.

The Darcy friction factor is defined as

where $τ12*$ is a component of the dimensional version of the elastic stress tensor τ and ⟨U^*⟩ is the cross-sectional average of U^*(r^*). For Newtonian flows, the same definition holds, but without the second term within the square brackets. The inequalities in Eq. (18) translate into the following inequalities for the friction-factor:

The right-most term of Eq. (20) is f for a Newtonian fluid with viscosity same as the total viscosity of the solution. The middle and the left-most terms are f for the polymer solution and the pure Newtonian solvent, respectively. These expressions for the friction-factor can be derived from Eqs. (6), (9), and (19) in the case of polymeric solutions and by using the parabolic profiles in the case of Newtonian flows.

To analyze the influence of Wi and L, we consider the limit β → 0 whereby the solutions are such that F = 1 + 8Wi²r²/L² and U(r) = 1 − r² + 4Wi²(1 − r⁴)/L². These relations show that the volume flow rate increases with Wi and decreases with L. Indeed, a large Wi implies that the relaxation time is far greater than the time scale of the flow, allowing the polymer strain—i.e., the stretching of the ends of the polymer molecules—to be increased by the mean shear. In what follows, we show that an enhanced stretching of polymers reduces the drag when compared with a Newtonian fluid with the same viscosity as the total viscosity of the non-Newtonian case. Moreover, a large L implies that the ratio of $Cii$ to L becomes small, resulting in a reduction of the restoring elastic modulus F. This, in turn, generates an increase in C₂₂ and C₃₃ given by Eq. (5). As shown below, an increase in C₂₂ yields an increase in the effective viscosity.

The profiles of the non-zero components of C_ij(r) are shown in Fig. 1(d). The large values of C₁₁ originate from the polymers undergoing stretching in the axial direction. On the other hand, the fact that C₂₂ and C₃₃ have values below one implies contractions in both the radial and azimuthal directions. The term C₁₂ plays a crucial role in translating polymer strain represented by C₁₁ into an enhancement of the strain rate U_r of the flow, as revealed by the ê₁ê₁ component of Eq. (3): 2C₁₂U_r = (FC₁₁ − 1)/Wi.

F(r) is essentially a representation of the elastic modulus [see Fig. 1(e)]. The polymers undergo maximum stretching near the wall given the higher values of the strain rate U_r. This causes an increase in F; a feature of FENE models that is absent in the Oldroyd-B model.

We can now explain the lower drag in the polymeric solution when compared to that of the Newtonian fluid with identical total viscosity. The divergence of the Newtonian stress tensor (⁠$∇$·τ) · ê₁ [see Fig. 1(f)] without polymer stands at a constant value of −4, with the negative sign implying (positive) dissipation. When polymers are added, a ratio of (1 − β) of this viscous dissipative component is replaced by the component due to polymer stretching, which is also dissipative of the fluid momentum. As evident from Fig. 1(f), this component is always greater than −4 throughout the flow field except at the center of the pipe. [Note that this component can be written as ${−2r2/[β(F−1)+1]}r/r$⁠, which takes the value of −4 in the limit r → 0, which is same as the Newtonian component.] To counterbalance the constant pressure gradient term, the Newtonian part readjusts itself to increase the dissipation by increasing the overall strain rate, thereby resulting in an enhanced volume flow rate. Such an increase in the strain rate is characterized by an increase in the slope of the velocity near wall with respect to 1 − r in Fig. 1(a).

To explain the decrease in |(⁠$∇$·τ) · ê₁| as r increases, we first note that $C12≡⟨R̃1R̃2⟩=0$ at the pipe center, thus implying that the polymer molecules undergo an uncorrelated random motion there. Since C₁₁ increases from the value of 1.0 with respect to r, the polymers undergo a rapid stretching in the x-direction [see Fig. 1(d)], causing the $R̃2$ component to pick-up a negative correlation with the $R̃1$ component that undergoes stretching. This negative correlation is due to the resistance to elongation associated with the restoring tendency F. However, this negative correlation starts to decrease with respect to r partly due to the slow-down in the axial stretching C₁₁ and the restoring tendency, F—i.e., C_11rr < 0 and F_rr < 0 for most part of the pipe. This causes less drag as the absolute value of the resistance term $[(C12F)r+C12F/r]/Wi$ comes down.

This phenomenon of C₁₁ or C₁₂ playing a crucial role in the laminar regime is in stark contrast with the turbulent regime, where the transverse diagonal component C₂₂ (or C₃₃) plays a critical role—it acquires a radial distribution similar to the effective viscosity.^18,24 In fact, both interpretations are reconciled as explained below.

The decrease in flow rate when compared to that of the pure solvent, or, the increase in the same when compared to that of a viscosity-matched Newtonian fluid, can also be understood by considering the effective viscosity μ_eff, which is defined as

Taking note of the facts that F → ∞ in the limits Wi → ∞ or L → 0 and F → 1 in the limits Wi → 0 or L → ∞, we arrive at the following bounds for μ_eff/μ:

which contains the analogous information as the inequalities in Eq. (18), i.e., the effective viscosity is less than the total viscosity, but higher than the viscosity of the solvent. In Fig. 2, μ_eff/μ is shown for a parameter set equal to that of Fig. 1(d). The figure confirms that the ratio μ_eff(r)/μ obeys the inequalities in Eq. (22). Clearly, μ_eff(r)/μ has a trend with respect to r, which is opposite to that of C₁₁. This shows that the stretching in the axial direction reduces the effective viscosity. As a general non-auxetic matter, the polymer exhibits an opposite trend in C₂₂ and C₃₃ by the way of contraction of the polymer dumbbells. This sets C₂₂ and C₃₃ to follow the same trend as μ_eff(r)/μ [see the inset of Fig. 1(d)]. As stated previously, this phenomenon is well-known in the turbulent regime.^18,24 The linear relation (up to the leading order) between μ_eff(r)/μ and C₂₂ is revealed by re-writing Eq. (21) as

where the second term within the braces is shown to be of O(r²).

In summary, the laminar regime of this non-Newtonian flow has interesting bounds for flow and friction determined by two corresponding Newtonian flows. The flow of a polymeric solution with a total viscosity identical to that of a Newtonian fluid exhibits lower drag on comparison with the latter. Nonetheless, the experienced drag is higher than that of the pure solvent under the same driving pressure gradient. It should be noted that these bounds are also valid for plane Poiseuille flow.

During the course of arriving at these bounds, we also solved the FENE-P model in a form widely known to the community researching turbulence and transition. These alternative solutions to those found by Cruz et al.²⁷ will enable the study of small perturbations to the FENE-P model expressed in Eqs. (1) and (2). However, any future study on transition should take into account that the stability of the viscoelastic flow at a certain Re should be compared against that of the corresponding flow of the solvent at Re/β in order to maintain the same pressure gradient. If not, the comparison would not be against the stability of the flow of the solvent, but with another Newtonian fluid with the same viscosity as the total viscosity of the solution.

This rich phenomenology can also be interpreted as an explanation of the well-known shear-thinning effect through the following sequence of comprehensive arguments: (1) The strain rate of the flow stretches the polymers in the axial direction, a phenomenon also known from the earlier literature; (2) such an extension causes contraction in the radial direction due to the resistance to elongation exerted by the restoring force between the ends of the polymer molecules; (3) this triggers a nonzero negative correlation between the x-component and r-component of the end-to-end polymer vector, resulting in an elastic shear stress (τ₁₂) to develop; (4) since there is no strain rate at the pipe center, this correlation is zero due to the liberty of the ends of the polymers to be at random motion; (5) due to the increase in the elastic modulus with respect to r, a feature of the FENE-P model, the stretching rate with respect to r slows down, giving rise to a convexity of C₁₁ (i.e., C_11rr < 0), which in turn causes concavity in τ₁₂; (6) The points in (3)–(5) imply that the divergence of elastic shear-stress, i.e., the resistive force contribution to the flow by the polymers decreases with r; (7) the driving force, i.e., the pressure gradient, which is balanced completely by the Newtonian viscous force and elastic force, is a constant; (8) since the elastic force decreases radially, the Newtonian viscous force should increase with r; (9) finally, this is achieved by enhancing the strain rate owing to the fact that the dynamic viscosity of the solvent is a constant, thus increasing the flow rate. The last point in this sequence could be perceived as shear-thinning by introducing a variable viscosity as a factoring function together with ∇²U, which balances the total dissipation.

Indeed, we proved that U and the components of the tensor C exhibit even or odd symmetries with respect to r. As an outlook, these symmetries can be exploited to derive conditions for the regularity of perturbations at the pipe center, as derived in the Newtonian case.³² These are expected to be valuable for the study of the stability and transition of viscoelastic pipe flows.

We thank Dr. Jurriaan J. J. Gillissen of the Department of Mathematics, University College London, for his fruitful comments on this letter.