Polymers in turbulence: any better than dumbbells?

Abstract

Polymer chains in turbulent flows are generally modelled as dumbbells, i.e. two beads joined by a nonlinear spring. The dumbbell only maps a single spatial configuration, described by the polymer end-to-end vector, thus a multi-bead FENE (finitely extensible nonlinear elastic) chain seems a natural improvement for a more accurate characterisation of the polymer spatial conformation. At a large Weissenberg number, a comparison with the more accurate Kuhn chain reveals that the multi-bead FENE chain drastically overestimates the probability of folded configurations. Surprisingly, the dumbbell turns out to be the only meaningful bead-spring model to coarse-grain a polymer macromolecule in turbulent pipe flows.

1. Introduction

Polymer dynamics in turbulent wall-bounded flows has gained interest in the last decades because of the capability of polymer chains to drastically modify the fluid dynamics of a Newtonian solvent. Polymer drag reduction (Benzi & Ching 2018), elastic turbulence (Steinberg 2021) and elasto-inertial turbulence (Dubief, Terrapon & Hof 2022) are three examples of the effect of the interaction between Newtonian solvent and polymers. Polymer drag reduction has been extensively investigated by means of direct numerical simulations (DNS) since the 1990s, with the polymer phase typically modelled as FENE-P (finitely extensible nonlinear elastic with Peterlin's approximation) dumbbells (Sureshkumar, Beris & Handler 1997; De Angelis, Casciola & Piva 2002), and only recently using Lagrangian FENE dumbbells to account for the polymer dynamics and effects in free space (Watanabe & Gotoh 2013) and wall-bounded (Serafini et al. 2022, 2023) flows. The dumbbell can only map the end-to-end distance of a polymer molecule, but none of the internal configurations that it can assume. Consequently, the polymer effect on the solvent can be only related to the polymer end-to-end distance. One of the open questions is whether or not the internal degrees of freedom of the polymer chain, e.g. spatial configuration, may have a role in the dynamics of the polymer suspension. In prototypal flows, e.g imposed uniaxial extension and laminar shear flows, the issue of coarse-graining from freely jointed chains to FENE dumbbells has been addressed (Doyle & Shaqfeh 1998; Somasi et al. 2002), while fewer studies discuss the adequacy of the FENE dumbbell model in turbulent flows. One-way coupled simulations of Zhou & Akhavan (2003) highlight that at a moderate Weissenberg number (ratio of the polymer relaxation time to the fluid time scale) an ad hoc FENE dumbbell can provide a reasonably good prediction of the a posteriori polymer stress in a turbulent wall-bounded flow when compared with a 10-bead FENE chain. Their results confirmed those obtained by Jin & Collins (2007), Watanabe & Gotoh (2010) and Vincenzi et al. (2021) in isotropic turbulence at comparable Weissenberg numbers. No investigations are available for large-Weissenberg-number polymer chains, which, however, characterise drag-reducing solutions of biological polymers.

In this paper, we use the multi-bead FENE model to investigate polymer dynamics in a turbulent pipe flow at a moderate friction Reynolds number $ {\textit {Re}}_{\tau }=180$ and large Weissenberg number $Wi=10^4$. The FENE model – both multi-bead and dumbbell – is compared with the Kuhn chain, which represents a thoroughly faithful model of an actual polymer macromolecule.

2. Methodology

The evolution of the Newtonian solvent in a pipe geometry is described by the (dimensionless) incompressible Navier–Stokes equations, completed with the impermeability and no-slip condition, $\boldsymbol {u}(t,\theta,r=1,z)=0$, at the wall:

(2.1a,b)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0,\quad \frac{\partial{\boldsymbol{u}}}{\partial{t}} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{u} \otimes \boldsymbol{u}) ={-}\boldsymbol{\nabla} p + \frac{1}{{\textit{Re}}} \boldsymbol{\nabla}^2 \boldsymbol{u}, \end{equation}

where $t$ is the time and $\boldsymbol {x}=\boldsymbol {x}(\theta,r,z)$ is the position vector function of the cylindrical coordinates, $\boldsymbol {u}(t,\boldsymbol {x})$ denotes the fluid velocity, $p(t,\boldsymbol {x})$ the hydrodynamic pressure, and $ {\textit {Re}}=U_{b}^*\ell _{0}^*/\nu ^*$ is the Reynolds number built with the pipe radius $\ell _0^*$, the bulk velocity $U_b^* = Q_b^*/({\rm \pi} \ell _0^{*2})$ ($Q_b^*$ is the volumetric flow rate), and the kinematic viscosity $\nu ^*$ (the asterisks denote dimensional quantities). The statistical steady-state flow rate is sustained by a constant pressure gradient in the axial direction, corresponding to a friction Reynolds number $ {\textit {Re}}_\tau =u_\tau ^* \ell _0^* /\nu ^*=180$ ($u_\tau ^*$ is the friction velocity). Periodic boundary conditions are applied in the axial and azimuthal directions.

Flexible polymers can take up a large number of configurations by the rotation of chemical bonds, therefore their shape can only be usefully described from a statistical point of view (Doi & Edwards 1988). To this purpose, a very simple model is the so-called freely jointed chain (FJC), which consists of $N$ beads connected by $N-1$ rodlike links of length $l$, each one able to rotate independently of the others. The Brownian dynamics of an FJC can be described by means of the overdamped Langevin equations (Cruz, Chinesta & Regnier 2012), which establishes the instantaneous mechanical equilibrium among the hydrodynamic, Brownian and the constraint forces needed to preserve the rod length $l$. The equation for the generic polymer chain reads

(2.2) \begin{align} \left.\begin{array}{l@{}} \displaystyle\frac{{\rm d}\kern0.7pt {\boldsymbol{x}_{c}}}{{\rm d}{t}} = \frac{1}{N}\sum_{n=1}^{N} \boldsymbol{u}_n + \frac{1}{N}\sum_{n=1}^{N} \boldsymbol{\xi}_{n} \\ \displaystyle \frac{{\rm d}{\boldsymbol{r}_{n}}}{{\rm d}{t}} = (\boldsymbol{u}_{n+1} - \boldsymbol{u}_{n}) +(f_{n-1}\boldsymbol{r}_{n-1}-2f_{n}\boldsymbol{r}_{n}+f_{n+1}\boldsymbol{r}_{n+1}) +\dfrac{r_{eq}(\boldsymbol{\xi}_{n+1}-\boldsymbol{\xi}_{n})}{\sqrt{3 Wi g(N)}} \end{array}\right\}, \end{align}

where the $n$th bead position is $\boldsymbol {x}_n$, $\boldsymbol {x}_{c}=\sum \boldsymbol {x}_{n}/N$ is the polymer centre, $\boldsymbol {r}_n=\boldsymbol {x}_{n+1}-\boldsymbol {x}_{n}$ is the link vector, and $\boldsymbol {u}_{n}$ is the fluid velocity at the position $\boldsymbol {x}_n$. The white noise terms $\boldsymbol {\xi }_{n}$ set the equilibrium polymer size $r_{eq}$ in a quiet solvent, $Wi=\tau ^*/(\ell _{0}^*/U_{b}^*)$ is the Weissenberg number of the polymer chain with $\tau ^*$ the polymer relaxation time, and $f_n$ are the constraint forces that preserve the length $l$ of each polymer segment. Given a polymer with contour length $L$ and equilibrium size $r_{eq}$, the FJC which has $N-1$ links of length $l = r_{eq}^2/L$ is called the Kuhn chain. The quantity $l_K=l$, called Kuhn length, accounts for the bending stiffness of the polymers on scales smaller than $l_K$ (Doi & Edwards 1988).

A further coarse-graining of a polymer consists of replacing the rods of the Kuhn chain with one, or more, entropic springs. The value of the entropic stiffness $k$ is that which sets the correct chain equilibrium size $r_{eq}$ in a quiet solvent. In the context of polymer dynamics, such chains are known as Gaussian chains. In particular, it can be proved that an FJC of $N$ beads (given $N$ large enough) can be coarse-grained into a chain of $\tilde {N}<< N$ beads where each single bond length has a Gaussian distribution (Doi & Edwards 1988). For a Gaussian chain, the local properties of a polymer are lost, but the global properties are well described on large scale. However, the chain, which can be also represented by a mechanical model that considers $N$ beads connected by a Hookean potential, presents the unrealistic feature that the polymer contour length $L$ can be overcome with a finite probability (Bird et al. 1987), and it is thus unsuitable for describing polymer dynamics in turbulent flows. This issue is typically avoided by replacing the Hookean spring law with a nonlinear law of empirical derivation, known as Wagner or FENE spring law (Bird et al. 1987). The FENE chain dynamics is still described by an overdamped Langevin equation (Picardo et al. 2020), which is formally the same as (2.2). In this case, $f_n$ are nonlinear functions of the link length, instead of the constraint forces. In particular, given the maximum link length $l$,

(2.3a,b)\begin{equation} f_n = \frac{1}{Wi_N} \left( 1- \frac{\| \boldsymbol{r}_n\|^2}{l^2} \right)^{{-}1},\quad Wi_N = \frac{2g(N)}{N-1}Wi. \end{equation}

In (2.2) and (2.3a,b), the function $g(N) = \gamma _N/\gamma _2$ is the ratio between the beads’ friction coefficient $\gamma _N$ of an $N$-bead chain (FJC or FENE) and the beads’ friction coefficient $\gamma _2$ of a FENE dumbbell with the same relaxation time. Figure 1 shows the trend of $g(N)$ for the FJC and the FENE chain. For the latter, the value of $g(N)$ can be directly evaluated, using (2.3a,b), once the link Weissenberg number $Wi_N$ is known. The value of $Wi_N$ is chosen such that the slowest eigenvalue of the (linearised) system (2.2) is $1/Wi$, as for the Rouse chain (Doi & Edwards 1988). For an FJC the friction coefficient $\gamma _N$ is computed by performing numerical simulations of the relaxation dynamics. The correct value of $\gamma _N$ is that for which the Weissenberg number of the chain results in being equal to the Weissenberg number $Wi$ of the polymer. In figure 1, we also show the trend of $g(N)$ derived from the Rouse scaling in the limit of large $N$, which for $N>10$ collapses on the FJC and FENE scaling. The system (2.1a,b) is numerically solved in cylindrical coordinates on a staggered grid by means of a second-order central scheme in space and a projection method that enforces the velocity divergence-free constraint. Time integration is performed using a four-step, third-order, low-storage Runge–Kutta scheme. The domain dimensionless size is $(2{\rm \pi} \times 1 \times 6{\rm \pi} )$, while the number of grid points are $[N_\theta \times N_r \times N_z] = [128 \times 129 \times 384]$. A coordinate transformation in the radial direction ensures at the wall a minimum grid size $\Delta r^+ \simeq 0.5$, and at the pipe centre a maximum grid size $\Delta r^+ \simeq 2$. The $+$ superscript denote dimensionless variables with respect to the wall unit $\ell _\tau ^*=\nu ^*/u_\tau ^*$.

Figure 1. The $g(N)$ for the FENE chain and FJC. Solid green line reports the analytic prediction given by the Rouse model, in case of $N\gg 1$.

With regards to the polymer system (2.2), time integration is performed with a two-step, second-order, low-storage Runge–Kutta scheme. When integrating the FJC, the constraint on the constant rod length is satisfied within a specified tolerance at each time step, see Liu (1989) for the detailed scheme. Polymer time step can be lower than that of the fluid, up to 100 times for the Kuhn chain with $N=201$, that is characterised by the fastest dynamics. Once the system reaches the statistical steady state, statistics are collected over the total number of polymers on a dimensionless time interval of $200$.

In all simulations, the polymer contour length is fixed to $L=0.04$ to mimic a $20\,\mathrm {\mu }{\rm m}$ double-stranded DNA macromolecule (Harnau & Reineker 1999). The Kuhn bond length is $l_K\approx 2\times 10^{-4}$ and the equilibrium length is $r_{eq}=2.83\times 10^{-3}$. For FENE chains, the equilibrium length is kept constant to $r_{eq}$ and the maximum link length is $l=L/(N-1)$. For FJCs, the bond length is again $l=L/(N-1)$. According to this definition, the equilibrium end-to-end distance of an FJC is equal to $r_{eq,FJC}=r_{eq}\sqrt {(N_k-1)/(N-1)}$ ($N_k$ is the number of beads of the Kuhn chain), since it is not possible to fix equilibrium and contour length at the same time. Since polymer chains are mostly fully extended in a turbulent flow at a large Weissenberg number, we believe that the physical relevant value to prescribe is the contour length $L$ rather than the equilibrium length in a quiet solvent.

A population of $N_p=10^5$ polymers with Weissenberg number $Wi=10^4$ is considered in the simulations. The chosen number of polymers is enough to get convergent statistics and allows the modification induced by polymers to the flow field to be safely disregarded. Despite at a large Weissenberg number, polymers can modify turbulence in a non-trivial way (Rosti, Perlekar & Mitra 2023); here the number of polymers is too small to induce any modification, being three orders of magnitude smaller than the amount needed to induce a turbulent modification, see Serafini et al. (2022) where $N_p=10^8$ polymers are employed to observe a significant backreaction.

3. Results

The results of the FENE dumbbell and FENE chain at $Wi=10^4$ are compared with a reference simulation of a Kuhn chain of $N=201$ beads, by reporting in figure 2(a) the probability density function (p.d.f.) of the end-to-end distance R normalised by L for the Kuhn chain (black line) and for the FENE chains with different numbers of beads (coloured lines). The p.d.f. provided by the FENE dumbbell model is very similar to that of the Kuhn chain, with a marked peak characterising the fully extended configuration. The FENE dumbbell slightly overestimates the probability of maximum stretching with respect to the reference Kuhn chain, for which the fully stretched state is characterised by all its internal links aligned. This configuration turns out to be the most probable and can be univocally mapped by the dumbbell end-to-end vector. All the intermediate values of the end-to-end distance can be instead associated with multiple configurations, and each one can in principle correspond to a different value of polymer stress, see Graham (2004). Thus it is natural to think about multiple-mode FENE chains to characterise internal configurations, without solving the entire dynamics of an FJC. Intuitively, increasing the number of beads, and thus the number of representable internal configurations, the p.d.f. given by the Kuhn chain should be progressively recovered. However, if $N\ge 3$, the extended state's probability drastically drops, and new intermediate peaks appear in the p.d.f.s. Since the majority of the internal links of the chains are fully stretched, see the p.d.f. of the normalised link length in figure 2(b), the new peaks in figure 2(a) can be associated with folded chains. For the 3-bead chain (red curve), the end-to-end distance p.d.f. has two peaks: one at $R/L=1$ and another at $R/L=0$ with fully extended internal links. A further increase in the number of beads results in the appearance of new peaks, with a spacing equal to $2/(N-1)$. Convergence is thus never reached by increasing the number of beads of a FENE chain. A more detailed description of the polymer shape can be provided by considering the joint p.d.f. of the end-to-end distance and the radius of gyration (Levi-Civita & Amaldi 1923)

(3.1)\begin{equation} R_g = \sqrt{\frac{1}{N} \sum_{n=1}^{N} \| \boldsymbol{x}_n-\boldsymbol{x}_c \|^2}, \end{equation}

a quantity that characterises the average volume occupied by a polymer chain. The joint p.d.f., for a 3-bead chain, is reported in figure 3(a) where both $R$ and $R_g$ are normalised by the polymer contour length $L$. All the possible configurations are bounded between the lines $R_g/L = R/(\sqrt {6}L)$ and $R_g/L = \sqrt {1/18+(R/L)^2/9}$. The former corresponds to two aligned links of equal length $r$, configuration A in figure 3(d), and the latter to two fully extended links of length $l$ misaligned by an angle $\alpha$, configuration B. Two possible configurations are characterised by zero end-to-end distance: $\boldsymbol {x}_1=\boldsymbol {x}_2=\boldsymbol {x}_3$ which has $R_g/L=0$, and $\boldsymbol {x}_1=\boldsymbol {x}_3\neq \boldsymbol {x}_2$, namely a chain with two superimposed links of equal length $\| \boldsymbol {r}_1 \|= \| \boldsymbol {r}_2 \| = r$, for which $R_g/L=\sqrt {2} r/(3L)$. In accordance with figure 2(a), the joint p.d.f. has two peaks, one at $R/L=0$ and $R_g=\sqrt {2}/6$ corresponding to the folded polymer (configuration B with $\alpha =0$), and one at $R/L=1$ and $R_g/L=1/\sqrt {6}$ corresponding to the fully extended polymer (configuration A with $r=l$). These two peaks are joined by a large probability line, denoted with the letter C in figure 3(a), that is well described by the curve

(3.2)\begin{equation} \frac{R_g}{L} = \sqrt{\frac{1}{9} \frac{R}{L} \left(\frac{R}{L}-1\right)+ \frac{1}{18}+\frac{1}{9}\frac{R^2}{L^2}}. \end{equation}

The curve has no direct analytic derivation, but is the simplest (quadratic) curve that connects the two peaks and is tangent to the line $R_g/L = R/(\sqrt {6}L)$ for $R/L=1$.

Figure 2. The p.d.f. of the normalised end-to-end distance $R/L$ (a), and of the link length $r/l$ (b), for FENE chains with different numbers of beads. The solid black line in (a) reports the p.d.f. for the Kuhn chain ($N=201$), taken as reference for all the other coarse-grained models.

Figure 3. (a–c) Joint p.d.f. of the normalised end-to-end distance $R/L$ and the normalised radius of gyration $R_g/L$ for the 3-bead FENE chain at $Wi=10^4,10,10^2$, respectively. (d) Sketch of two reference polymer configurations. Configuration A represents a 3-bead chain with two aligned links with same length $r\le l$. Configuration B represents a 3-bead chain with fully extended links misaligned by an angle $\alpha$.

The weird behaviour of the FENE chain is peculiar to polymers with high Weissenberg number. The joint p.d.f. of end-to-end distance and radius of gyration is reported at $Wi=10$ and $Wi=10^2$ in figure 3(b,c), respectively. At $Wi=10$ the probability is concentrated very close to the line $R_g/L = R/(\sqrt {6}L)$, i.e. configuration A in figure 3(b), while at $Wi=10^2$ the p.d.f. is qualitatively similar to $Wi=10^4$, with a less marked peak for the folded configuration. A possible interpretation of the difference between the two cases is the following. At moderate $Wi$ the entropic relaxation is not negligible in determining the polymer extension, thus, when the turbulent flow field is not stretching the polymer, the internal links tend to relax remaining aligned. At large $Wi$, instead the polymer is a ‘slave’ of the turbulent flow field, which controls both the macromolecule's extension and relaxation dynamics. The elastic restoring forces are always negligible with respect to the hydrodynamic contribution and are significant only in the fully extended configuration to avoid polymers overcoming their maximum length. When the turbulent flow field is not extending the polymer, the internal links tend to misalign themselves, sampling with large probability configurations on the curve C in figure 3(a) that leads to the fully extended folded state.

The joint p.d.f. of end-to-end distance and radius of gyration is also reported in figure 4 for the Kuhn chain (a), and the 21-bead FENE chain (b). For the Kuhn chain, the configurations with $R/L=0$ are never found, and the maximum probability is associated with the fully elongated polymer. Besides a non-negligible probability of the $R/L=0$ configurations, the 21-bead chain also samples folded states with fully extended links, with a constant spacing determined by the maximum link length $2/(N-1)$, see figure 4(b). At a large Weissenberg number, any increase in the number of beads in a FENE chain results in a worse description of the polymer geometry, as the chain samples a larger number of configurations with multiple ways folded and fully extended links. Surprisingly, the most accurate FENE chain is the dumbbell that, being only able to describe the polymer extension, cannot sample any of the improbable states that the multiple-mode chains describe. In the plane $(R/L, R_g/L)$, the dumbbell can only sample configuration on the line passing through the origin, corresponding to configuration A in figure 3(b). Remarkably, this is the region that the Kuhn chain samples with the largest probability.

Figure 4. Joint p.d.f. of the normalised end-to-end distance $R/L$ and the normalised radius of gyration $R_g/L$ for the Kuhn chain ($N=201$) (a), and for the FENE chain with $N=21$ (b), at Weissenberg number $Wi=10^4$.

Based on the above results, the only meaningful model to fine-grain the dumbbell in turbulent pipe flows is the Kuhn chain, which, however, requires an enormous computational effort since 201 beads are needed to reproduce the actual polymer properties at equilibrium. In contrast, in a turbulent flow, polymers are likely to be very far from the equilibrium condition, being most of the time stretched by the turbulent fluctuations. Since the Kuhn bond length is approximately 30 times smaller than the smallest turbulent scale, we can conjecture that an FJC with fewer beads, thus a different (larger) equilibrium size, can still adequately describe the polymer dynamics in a turbulent flow, when the Weissenberg number is large. Figure 5 shows the p.d.f. of the end-to-end distance of FJCs with different numbers of beads. The p.d.f.s are almost indistinguishable if $N\ge 21$, thus a 21-bead FJC can effectively characterise the polymer dynamics of the considered DNA macromolecule at $ {\textit {Re}}_\tau =180$ and $Wi=10^4$. Remarkably, the bond length of a 21-bead chain is approximately half of the wall unit. With the number of beads of a 21-bead chain being 10 times smaller than the number of beads of the Kuhn chain, the time required for the solution of system (2.2) can be reduced by a factor of $10$ at least. For a smaller number of beads, the FJC samples some of the folded configurations, as in the case of the multi-bead FENE chain. For a 3-bead chain, the folded configuration appears with a non-negligible probability and in the plane $(R/L, R_g/L)$ only configuration B can be sampled. During its relaxation, a 3-bead FJC cannot move on line C and the probability of zero end-to-end distance is smaller than in the case of a 3-bead spring chain. In conclusion for an FJC, increasing the number of links (thus the number of internal configurations), folded states become progressively less probable and convergence is reached.

Figure 5. The p.d.f. of the normalised end-to-end distance $R/L$ for FJCs with different numbers of beads at Weissenberg number $Wi=10^4$.

4. Conclusions

In this paper, we have compared the FENE dumbbell model with the multi-bead FENE and freely-jointed chain at $Wi=10^4$ in a turbulent pipe flow.

The comparison of the coarser-grained model with a reference Kuhn chain reveals that the multi-bead FENE chain is an improper characterisation of the polymer dynamics in a turbulent pipe flow at large Weissenberg number and convergence is never reached by increasing the number of internal degrees of freedom. Surprisingly, the only effective model to coarse-grain a Kuhn chain is the FENE dumbbell, while multi-bead FENE chains sample spatial configurations with folded fully extended links. These folded configurations are demonstrated to be improbable by direct comparison with the Kuhn chain. Conversely, for FJCs, the increase in the number of beads ensures convergence to the statistics of the Kuhn chain. With the entropic relaxation being negligible at large Weissnberg number, in the present condition, the polymer dynamics of the DNA molecule is well captured by a 21-bead FJC with a speed-up of order 10 in the computational time. We remark that a 21-bead FJC no longer describes the large-scale properties of a polymer in a quiet solvent, but it is sufficiently accurate to capture the polymer dynamics in a turbulent flow, with the polymer chain always being far from the equilibrium.

In conclusion, the dumbbell model is accurate in the description of the polymer extension statistics in turbulent pipe flows. Nonetheless, if the internal degrees of freedom of the chain must be accounted for, the only meaningful model is the FJC with bond length at most comparable to the smallest relevant hydrodynamical scale. We expect these conclusions to apply to shear-dominated flows, while other kinds of dynamics, such as those occurring in free-space conditions, deserve further investigations at Weissenberg numbers comparable to the present study.