Signatures of elastoviscous buckling in the dilute rheology of stiff polymers
Brato Chakrabarti, Yanan Liu, Olivia du Roure, Anke Lindner, David Saintillan
SSignatures of elastoviscous buckling in the dilute rheology of stiff polymers
Brato Chakrabarti,
1, 2
Yanan Liu, ∗ Olivia du Roure, Anke Lindner, and David Saintillan † Department of Mechanical and Aerospace Engineering,University of California San Diego, 9500 Gilman Drive, La Jolla, CA 92093, USA Center for Computational Biology, Flatiron Institute, New York, NY 10010, USA Laboratoire de Physique et M´ecanique des Milieux H´et`erog`enes, UMR 7636,ESPCI Paris, PSL Research University, CNRS, Universit´e de Paris,Sorbonne Universit´e, 10 rue Vauquelin, 75005 Paris, France
As a stiff polymer tumbles in shear flow, it experiences compressive viscous forces that can causeit to buckle and undergo a sequence of morphological transitions with increasing flow strength. Weuse numerical simulations to uncover the effects of these transitions on the steady shear rheology ofa dilute suspension of stiff polymers. Our results agree with classic scalings for Brownian rods inrelatively weak flows but depart from them above the buckling threshold. Signatures of elastoviscousbuckling include enhanced shear thinning and an increase in the magnitude of normal stress differ-ences. We discuss our findings in the light of past work on rigid Brownian rods and non-Brownianelastic fibers and highlight the subtle role of thermal fluctuations in triggering instabilities.
I. INTRODUCTION
Understanding and relating bulk rheological properties of complex fluids to the orientations, deformations, andinteractions of their microscopic constituents has been a long-standing challenge in fluid mechanics [18, 23]. With theability to visualize single molecules using fluorescence microscopy, a large body of research using deoxyribonucleicacid (DNA) has focused over the last two decades on deciphering the dynamics and rheological properties of dilutelong-chain polymer solutions in simple flows [26]. The persistence length (cid:96) p of DNA molecules is much shorter thantheir typical contour length L , and in this limit conformational properties are governed by a competition betweenentropic forces favoring coiled conformations and viscous forces tending to stretch the molecules, as quantified by theWeissenberg number W i ≡ ˙ γτ r , or product of the applied strain rate ˙ γ with the longest polymer relaxation time τ r .Along with coarse-grained simulations [15] and kinetic theories [28], these molecular rheology studies have highlightedhow microscopic conformational dynamics give rise to macroscopic rheological properties such as shear thinning andnormal stress differences in shear flow [12, 25], and lead to the ‘coil-stretch’ transition in extensional flow [24].In the other limit, the rheology of rigid Brownian rod-like suspensions is also well understood [6]. Thermal fluctu-ations in this case result in orientational diffusion with rotational diffusivity d r = 3 k B T ln(2 r ) /πµL , where k B T isthe thermal energy, µ is the shear viscosity of the solvent, and r = L/a is the aspect ratio of the rods with character-istic radius a . The competition between orientational diffusion favoring an isotropic distribution and the backgroundshear that tries to align the rods is characterized by the rotary P´eclet number P e ≡ ˙ γd − r , where d − r describes theorientational relaxation time. With increasing P e , the preferential alignment of these rod-like polymers near the flowaxis reduces viscous dissipation, resulting in shear thinning and also yielding nonzero normal stress differences. Threedistinct scaling regimes for the shear viscosity and normal stress differences as functions of
P e have been identifiedand characterized by the foundational theoretical analyses of Leal & Hinch [19], Hinch & Leal [10] and Brenner [3](see § II C and Table I for a summary). However, so far only rigid Brownian rods have been considered in detail, andthe role of flow-induced deformations on the rheology of these suspensions remains largely unexplored.We address this problem here with focus on stiff polymers characterized by (cid:96) p (cid:29) L , the opposite limit comparedto DNA. While in weak flows these filaments behave as rigid rods, they are known to undergo various bucklinginstabilities in stronger flows [2, 4, 8, 20, 21], yet clear insight into the specific role of these instabilities in the rheologyof dilute suspensions is lacking. Here, we use numerical simulations to relate the morphological transitions of stiffBrownian filaments in simple shear flow to the rheology of their dilute suspensions. We also contrast our predictions ∗ Present address: School of Physics, Northwest University, Xi’an, China † Correspondence: [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] F e b with known results for non-Brownian deformable fibers [2, 27] and uncover how they are altered by shape fluctuationsand orientational diffusion.The paper is organized as follows. In § II, we provide details of the polymer model, measures of the extra stressand a brief summary of the scaling laws for Brownian rigid rods. We present numerical results for the rheology inboth two and three dimensions in § III and discuss our predictions in the context of the rheology of rigid Brownianrods as well as non-Brownian elastic fibers.
II. PROBLEM DESCRIPTION AND METHODOLOGYA. Governing equations
In the dilute limit, we simulate the dynamics of a single polymer modeled as a fluctuating, inextensible Euler elasticawith centerline parameterized as x ( s , t ) where s is arc-length [20]. Hydrodynamics is captured by local slender bodytheory for Stokes flow, in which the centerline position evolves as8 πµ [ ∂ t x ( s , t ) − u ∞ ] = − Λ · f ( s , t ). (1)Here, u ∞ = ( ˙ γy , 0, 0) is the background shear flow with constant shear rate ˙ γ . The force per unit length exerted bythe filament on the fluid is modeled as f = B x ssss − ( σ x s ) s + f b , where B is the bending rigidity, σ is a Lagrangemultiplier that enforces inextensibility of the filament and can be interpreted as line tension, and f b is the Brownianforce density obeying the fluctuation-dissipation theorem. The local mobility operator Λ accounts for drag anisotropyand is given by Λ · f = [(2 − c ) I − ( c + 2) x s x s ] · f , (2)where c = ln( (cid:15) e) < (cid:15) = r − is the inverse aspect ratio. We scale lengthsby L , time by the characteristic relaxation time of bending deformations τ r = 8 πµL /B , elastic forces by the bendingforce scale B/L , and Brownian forces by √ L/(cid:96) p B/L . The dimensionless equation of motion then reads ∂ t x ( s , t ) = ¯ µ u ∞ − Λ · (cid:104) x ssss − ( σ x s ) s + (cid:112) L/(cid:96) p ζ (cid:105) , (3)where ζ is a Gaussian random vector with zero mean and unit variance. Two dimensionless groups appear: (i) theelastoviscous number ¯ µ ≡ ˙ γτ r c = 8 πµ ˙ γL Bc (4)serves as the measure of hydrodynamic forcing against internal elasticity and plays a role analogous to the Weissenbergnumber for flexible polymers, and (ii) L/(cid:96) p captures the importance of thermal shape fluctuations. The limit of rigidBrownian rods formally corresponds to ¯ µ → L/(cid:96) p → µ and L/(cid:96) p are related to the rotary P´eclet number: P e = ¯ µ c
24 ln(2 r ) (cid:96) p L , (5)which will facilitate comparisons of our simulations of Brownian polymers with finite bending resistance with analyticalpredictions for rigid Brownian rods. B. Measures of stress
A calculation of the extra stress in a dilute suspension of force- and torque-free particles was provided by Batchelor[1]. The single-particle contribution to the bulk stress tensor in the dilute limit is given by the stresslet, whichgeneralizes the Kirkwood formula commonly used for molecular systems [13]. For our polymer model, the expressionfor the extra stress is Σ = − n (cid:68) (cid:90) L (cid:20)
12 ( xf + f x ) − I ( x · f ) (cid:21) d s (cid:69) , (6)where n is the number density in the suspension, f is the dimensionless force density exerted on the fluid withcontributions from both elastic deformations and Brownian fluctuations, and (cid:104)·(cid:105) denotes the ensemble average. This Regimes Scalings η Σ xy Σ xx , Σ yy , Σ zz P e (cid:28)
P e P e r + r − (cid:29) P e (cid:29) P e − / P e / P e / P e (cid:29) r + r − constant P e P e TABLE I. Asymptotic scalings for the relative viscosity η and the individual stress components in dilute suspensions of rigidBrownian rods in various regimes of the rotational P´eclet number P e and aspect ratio r [3, 10]. The first and second normalstress differences N > N < expression is extremely convenient for non-Brownian fibers [2]. However, in simulations of Brownian polymers,fluctuations have contributions of O (∆ t − / ) where ∆ t is the integration time step. These contributions also enterthe Lagrange multiplier σ ( s , t ) that enforces inextensibility, resulting in a poor convergence of the ensemble average asfirst pointed out by Doyle et al. [7]. Since our interest is in the steady-state extra stress, we instead use the Giesekusstress expression commonly used for polymers [7, 22], Σ = − n (cid:68) (cid:90) L (cid:20)
12 ( x R · u ∞ + R · u ∞ x ) − I ( x · R · u ∞ ) (cid:21) d s (cid:69) , (7)where R = Λ − is the local resistance tensor along the centerline. Results obtained with equation (7) were testedagainst (6) in various regimes of P e . In weak flows (
P e (cid:28)
P e were obtained by both methods across all regimes of
P e . In the following,we present results based on equation (7), which is computationally more efficient than the Kirkwood expression.
C. Summary of rigid rod rheology
A slender non-Brownian rod in shear flow undergoes a periodic tumbling motion known as a Jeffery orbit [14].During this periodic tumbling, the particle spends most of its time aligned with the flow direction and equal amountsof time in the extensional and compressional quadrants of the flow. This dynamics is fundamentally altered in thepresence of rotational diffusion. While the shear flow still results in quasi-periodic tumbling, the Brownian rod is nowable to stochastically sample different orbits [29], resulting in an anisotropic orientational probability distribution ψ ( p ) at steady state, where p is a unit vector that identifies the orientation of the rod [5]. This distribution leadsto a mean orientation of the rod in the extensional quadrant, giving rise to a contractile stresslet as the inextensiblerod resists stretching by the flow. This stresslet in turn alters the effective viscosity of the system. In the dilutelimit of nL (cid:28)
1, computing the extra stress reduces to obtaining the steady state orientation distribution of a singleBrownian rod. Contributions to the stresslet arise from the external flow and from Brownian diffusion and can becomputed using slender body theory [1, 3, 10, 19]: Σ f = πµnL r ) (cid:20) (cid:104) pppp (cid:105) − I (cid:104) pp (cid:105) (cid:21) : E ∞ , Σ b = 3 nk B T (cid:20) (cid:104) pp (cid:105) − I (cid:21) , (8)where E ∞ is the rate-of-strain tensor of the applied flow u ∞ . The extra stress Σ = Σ f + Σ b is thus entirely determinedfrom the second and fourth moments, (cid:104) pp (cid:105) and (cid:104) pppp (cid:105) , of rod orientations. These moments can be computed fromthe steady state orientation distribution function ψ ( p ), which is set by the balance of the advective rotational flux dueto the flow and of the Brownian diffusive flux. Hinch & Leal [10] solved for the distribution function and associatedparticle stress in three distinct asymptotic regimes of the rotary P´eclet number P e . These regimes and correspondingscalings are summarized in Table I. The three rheological measures of primary interest to us are the relative polymerviscosity η and the first and second normal stress differences N and N , η = Σ xy µ ˙ γnL , N = Σ xx − Σ yy nk B T , N = Σ yy − Σ zz nk B T . (9) -3 -2 -1 -2 -1 −
P e (cid:28)
1) and green (1 (cid:28)
P e (cid:28) r + r − ) regions. The vertical dashed lines show the onsets of C buckling and U turns [20]. ( f ) Typical sequences of conformations during tumbling, C buckling, and a U turn. Corresponding regimes andvalues of P e are labeled in ( a ). In all figures, marked scalings before the buckling transition are rigid rod predictions (solidline), whereas scalings past the transition are numerical observations (dotted line). III. NUMERICAL RESULTS AND DISCUSSIONA. Three-dimensional rheology of stiff polymers
We present numerical results on the rheology in three dimensions, with focus on the case of stiff slender polymerswith (cid:96) p /L = 1000 and aspect ratio r = 220. In this limit, the main effect of Brownian motion is to cause orientationaldiffusion with negligible shape fluctuations, and any deformations are thus the result of elastoviscous buckling. Fig-ure 1( a , b , c ) shows the relative polymer viscosity η and first and second normal stress differences N as functionsof P´eclet number P e (or, equivalently, elastoviscous number ¯ µ ). In weak flows ( P e (cid:28)
1, pink region), η exihibitsa plateau whereas N both grow from zero as P e , with N > N <
0. With increasing flow strength,shear-thinning takes place as the polymers start aligning with the flow, and a second regime begins with scalings of η ∼ P e − / and N ∼ P e / in perfect agreement with the theoretical predictions of Table I. The data for the secondnormal stress difference is very noisy in this range of P e and fails to capture the expected scaling for reasons weexplain below. In this regime, the filament remains straight and tumbles quasi-periodically as shown in the first rowof figure 1( f ).As the filament performs a tumble, it rotates across the compressional quadrant of the flow, where it experiencescompressive viscous stresses. Above a critical value of the elastoviscous number of ¯ µ (1) = 306.6, these stresses canovercome bending rigidity and drive an Euler buckling instability leading to deformed configurations reminiscent ofa C shape, typical of the first mode of buckling [2]. The filament then rotates as a C and stretches out again as itenters and sweeps through the extensional quadrant. With increasing flow strength, the filament becomes more likelyto buckle while remaining nearly aligned with the flow direction, and this ultimately gives rise to distinctive folded U shaped conformations that perform tank-treading motions while maintaining a mean orientation close to the flowaxis [9]. As uncovered in our past work [20], the transition to this new mode of transport occurs at ¯ µ (2) = 1.8 × .Both of these thresholds are indicated by vertical lines in figure 1( a , b , c ) and, for the chosen values of (cid:96) p /L and r , fallwithin the intermediate scaling range of 1 (cid:28) P e (cid:28) r + r − (green region). Typical conformations from tumbling, C buckling and U turns are shown in figure 1( f ).Quite remarkably, we find that the onset of buckling has no immediate signature on the rheology, with the Brownianrod-like scaling laws persisting past the critical value of ¯ µ (1) . This result is in contrast with the two-dimensionalrheology of non-Brownian elastic fibers studied by [2] and [27], where buckling is responsible for shear thinning andnonzero normal stress differences. This discrepancy is attributed to the presence of 3D rotational diffusion in oursimulations. In shear flow, the viscous compressive force experienced by the filament is a function of its orientation andreaches a maximum at an angle of π/ ψ ( p ), and therefore the probabilityof a buckling event is negligible at ¯ µ = ¯ µ (1) . As ¯ µ is increased beyond ¯ µ (1) , buckling becomes increasingly more likely,and indeed η and N start to depart from the intermediate scalings before the onset of tank-treading in figure 1.This departure marks a transition to new scalings of η ∼ P e − ± and N ∼ P e ± and is accompanied by asharp increase in N . These rheological changes are clear signatures of elastoviscous buckling, as they occur withinthe range of validity of the intermediate rigid rod scalings.A more complete picture is provided in figure 1( d , e ), showing the shear and diagonal components of the extra stresstensor. In particular, we find that normal stresses in figure 1( e ) are dominated by Σ xx , while Σ yy and Σ zz , which arenegative, have smaller magnitudes. All three components follows the same scaling, with the intermediate scaling of P e / giving way to a nearly linear scaling into the buckling regime. Note that for P e (cid:29) yy and Σ zz are almost identical, which explains the strong noise in the data for N in figure 1( c ), especially in the intermediatescaling regime.To relate these findings to filament conformations, we introduce the gyration tensor G ( t ) = (cid:90) ( x − x c ) ( x − x c ) d s , (10)whose dominant eigenvector is used to define the mean filament orientation θ ( t ) with respect to the flow direction. Itsensemble average is shown as a function of P´eclet number in figure 2( a ). In weak flows ( P e (cid:28) (cid:104) θ (cid:105) asymptotes tothe value of π/ (cid:104) G (cid:105) , whichdescribe the variance of the polymer mass distribution along the coordinate directions, all tend to the same value of1 / π in figure 2( b ). Alignment of the filament with the flow is accompanied by a decrease in the mean orientationangle (cid:104) θ (cid:105) with increasing shear rate, and is also indicated by the growth and saturation of (cid:104) G xx (cid:105) while (cid:104) G yy (cid:105) rapidlydecays. Increasing flow strength also forces the filament towards the shear plane leading to a decrease in (cid:104) G zz (cid:105) aswell. In strong flows, the initiation of U turns leads to a sharper decrease in both (cid:104) G yy (cid:105) and (cid:104) G zz (cid:105) , as the emergentfolded conformations remain increasingly aligned with the flow direction as seen in the third row of figure 1( f ).It is interesting to note that the scaling law for the viscosity is altered at a slightly lower value of ¯ µ compared to thegyration tensor or mean orientation angle. This solidifies the idea that the signature of deformations observed in the B u c k li n g U - T u r n -2 -1 −
Above, we have discussed the case of stiff polymers and have compared our results to known scaling laws for rigidrods in three dimensions. Stiff polymers and rods both experience strong rotational diffusion. Shape fluctuationsare absent in the case of rigid rods and remain weak for the stiff polymers in our simulations performed in the limitof (cid:96) p /L (cid:29)
1. The main difference between the two systems is thus the occurrence of buckling instabilities above agiven threshold for the stiff polymers. Numerical results show a clear signature of such elastoviscous buckling on therheology, with enhanced shear thinning and normal stress differences compared to the case of rigid rods.We now address the role of shape fluctuations on the rheology. Remaining in the limit of weak fluctuations forstiff polymers, we can tune fluctuations by varying (cid:96) p /L while keeping rotational diffusion important, with a smaller (cid:96) p /L corresponding to stronger shape fluctuations. For simplicity, we present here results obtained in 2D, and showin figure 3 the variation of the shear viscosity η and normal stress difference N as functions of shear rate for twocombinations of ( (cid:96) p /L , r ). Scaling laws have been determined for both η and N before and after the bucklingthreshold, with a transition region where the scaling is changing continuously. Before the threshold, the measuredscalings are in agreement with 3D rigid rod predictions within error bars, and quantitatively identical results areobtained for both values of (cid:96) p /L since elasticity plays a negligible role in that regime. As the P´eclet number increases,buckling occurs first for the smaller value of (cid:96) p /L , since P e and (cid:96) p /L are related by equation (5) and the bucklingthreshold occurs at a fixed ¯ µ = ¯ µ (1) . Beyond the threshold, close agreement is found between 2D and 3D results,and identical scalings are obtained for η and N for both values of (cid:96) p /L , indicating that varying the importance ofshape fluctuations, while remaining in the limit of weak fluctuations, does not alter the observed rheology noticeably.The limit of strong shape fluctuations, relevant to a number of experimental systems and to previous work [9, 20], isoutside the theoretical and numerical framework developed here.Another observation can be made from figure 3. In three dimensions in figure 1, the changes in the scaling lawsof viscosity and normal stress differences do not occur right at the theoretical onset of buckling, but are delayeddue to the presence of rotational diffusion, as the probability for the polymer to align perfectly with the directionof maximum compression is negligible at the buckling threshold. On the contrary, in 2D, any polymer performing atumble is required to sweep through the entire compressional quadrant and is therefore more likely to buckle. Thisinterpretation is confirmed by our 2D resuits in figure 3, in which the change in scaling due to deformations now occursslightly closer to the theoretical buckling threshold. Nevertheless, and in contrast to fully non-Brownian simulations[2], the transition is not abrupt, likely due to the presence of Brownian noise. The fact that the slope changes occurcloser to the buckling threshold confirms again that the occurrence of C buckling, rather than U turns, is responsible −
P e
The authors thank M. Shelley for useful discussions. A.L., B.C. and Y.L. acknowledge funding from the ERCConsolidator Grant PaDyFlow (Agreement 682367). D.S. acknowledges funding from a Paris Sciences Chair fromESPCI Paris. [1]
Batchelor, G. K.
J. Fluid Mech. , 545–570.[2] Becker, L. E. & Shelley, M. J.
Phys. Rev. Lett. , 198301.[3] Brenner, H.
Intl. J. Multiphase Flow ,195–341.[4] Chakrabarti, B., Liu, Y., LaGrone, J., Cortez, R., Fauci, L., du Roure, O., Saintillan, D. & Lindner, A.
Nat. Phys. , 689–694.[5] Chen, S. B. & Jiang, L.
Phys.Fluids , 2878–2890.[6] Doi, M. & Edwards, S. F.
The Theory of Polymer Dynamics . Oxford University Press.[7]
Doyle, P. S., Shaqfeh, E. S. G. & Gast, A. P.
J. Fluid Mech. , 251–291.[8] du Roure, O., Lindner, A., Nazockdast, E. & Shelley, M.
Annu. Rev. Fluid Mech. , 539–572.[9] Harasim, M., Wunderlich, B., Peleg, O., Kr¨oger, M. & Bausch, A. R.
Phys. Rev. Lett. , 108302.[10]
Hinch, E. J. & Leal, L. G.
J. Fluid Mech. , 683–712.[11] Huber, B., Harasim, M., Wunderlich, B., Kr¨oger, M. & Bausch, A. R.
ACS Macro Lett. , 136–140.[12] Hur, J. S., Shaqfeh, E. S. G. & Larson, R. G.
J. Rheol. , 713–742.[13] Irving, J. H. & Kirkwood, J. G.
J. Chem. Phys. , 817–829.[14] Jeffery, G. B.
Proc. R. Soc. London A , 161–179.[15]
Jendrejack, R. M., de Pablo, J. J. & Graham, M. D.
J. Chem. Phys. , 7752–7759. [16]
Kirchenbuechler, I., Guu, D., Kuraniawan, N. A., Koenderink, G. H. & Lettinga, M. P.
Nat. Commun. , 5060.[17] Lang, C., Kohlbrecher, J., Porcar, L., Radulescu, A., Sellinghoff, K., Dhont, J. K. G. & Lettinga, M. P.
Macromolecules , 9604–9612.[18] Larson, R. G.
J. Rheol. , 1–70.[19] Leal, L. G. & Hinch, E. J.
J. Fluid Mech. ,685–703.[20] Liu, Y., Chakrabarti, B., Saintillan, D., Lindner, A. & du Roure, O.
Proc. Natl. Acad. Sci. USA , 9438–9443.[21]
Manikantan, H. & Saintillan, D.
Phys. Rev. E , 041002.[22] ¨Ottinger, H. C. Stochastic Processes in PolymericFluids , pp. 1–15. Springer.[23]
Schroeder, C. M.
J. Rheol. , 371–403.[24] Schroeder, C. M., Babcock, H. P., Shaqfeh, E. S. G. & Chu, S.
Science , 1515–1519.[25]
Schroeder, C. M., Teixeira, R. E., Shaqfeh, E. S. G. & Chu, S.
Macromolecules , 1967–1978.[26] Shaqfeh, E. S. G.
J. Non-Newtonian Fluid Mech. , 1–28.[27]
Tornberg, A.-K. & Shelley, M. J.
J.Comput. Phys. , 8–40.[28]
Winkler, R. G.
Phys. Rev. Lett. , 128301.[29] Z¨ottl, A., Klop, K. E., Balin, A. K., Gao, Y., Yeomans, J. M. & Aarts, D. G. A.L.
Soft Matter15