Dynamics of dissolved polymer chains in isotropic turbulence

Polymers are remarkable molecules that have relaxation times that can span 15 orders of magnitude. The very longest of the relaxation times for high molecular weight polymers are sufficiently long to overlap with fluid mechanical times scales; under those circumstances, polymers can influence the flow. A well-known example that is still not fully understood is polymer drag reduction. It has been known since Toms (1949 Proc. 1st Int. Congress on Rheology 2 135–41) that parts per million (mass basis) concentrations of polymers can reduce the drag on a solid surface by as much as 80%. Understanding the mechanism of drag reduction requires an understanding of the dynamics of the dissolved polymer chain in response to local fluctuations in the turbulent flow field. We investigate this by using Brownian dynamics simulations of bead-spring models of polymers immersed in a turbulent solvent that is separately computed using direct numerical simulations. We observe that polymer chains with parameters that are effective for drag reduction generally remain stretched for long periods of time and only occasionally relax. The relatively restricted configuration space they sample makes it reasonable to represent their behavior with simpler dumbbell models. We also study the spatial structure of the polymer stresses using a Lagrangian strategy. The results explain the need for relatively high spatial resolution for numerical simulations of polymer flows.


Introduction
Polymers are generally long-chain hydrocarbons that are made up of a sequence of thousands or millions of smaller monomer units. The functionality and properties of the polymer are controlled by the physical/chemical properties of the monomer units, as well as the length and configuration (e.g. the degree of branching) of the chain. While the timescale of each monomer unit may be on the order of molecular times (e.g. 10 −15 -10 −12 s), the relaxation time for a high molecular weight polymer molecule can be on the order of milliseconds to seconds. The overlap of the longest relaxation time of the molecule with hydrodynamic time scales gives rise to viscoelastic and other non-Newtonian phenomena. One example of this is turbulent drag reduction in polymer solutions. It has been known since the early work of Toms (1949) that even parts per million (on a mass basis) concentrations of high-molecular-weight polymers in a Newtonian solvent can reduce the drag on a surface by as much as 80%.
There have been some important advances in our fundamental understanding of polymer-turbulence interactions that lead to drag reduction (L'vov et al 2004(L'vov et al , Vincenzi et al 2007. This has, in part, been driven by recent direct numerical simulations (DNS) of polymer drag reduction in channel flows (Dimitropoulos et al 1998, 2001, Ilg et al 2002, Min et al 2003a, 2003b, Ptasinski et al 2003, Sureshkumar et al 1997, boundary layers (Dimitropoulos 2005(Dimitropoulos , 2006 and homogeneous turbulent flows (De Angelis et al 2005, Vaithianathan and Collins 2003, Vaithianathan et al 2006, 2007. Numerical simulations have the advantage of simultaneously solving for the flow and polymer fields, allowing for a more detailed investigation of the mechanisms underlying the drag reduction phenomenon. The computational demand of three-dimensional DNS of turbulence limits the models for the polymer to a class known as the finitely extensible nonlinear elastic model that invokes the Peterlin (1961) approximation to close the dynamic equations (FENE-P). In this description, the polymer is modeled as two beads connected by a spring with a nonlinear force law that diverges as the separation distance approaches the maximum extension, L max (Beris and Edwards 1994, Bird et al 1987, Larson et al 1999. The model reproduces shear thinning, and yields an extensional viscosity that is substantially larger than the shear viscosity, which appears to be an essential ingredient to capture the qualitative features of turbulent drag reduction in wall-bounded flows. Dumbbell models such as FENE-P are only capable of representing the longest of the relaxation times of the polymer. However, as noted above, real polymers possess a broad spectrum of time scales associated with the relaxation of different segments of the molecule. This shortcoming of the FENE-P model can be partially addressed through multi-bead finitely extensible nonlinear elastic (FENE) formulations. Here, the dumbbell is replaced by a sequence of beads connected by nonlinear springs. The different relaxation rates associated with different subsets of the chain give rise to a spectrum of relaxation times that more closely mimics a real polymer molecule. These models are generally more accurate than dumbbell models; however, the computational expense is substantially higher, making them prohibitive in all but the simplest of flows (Laso andOttinger 1993, Ottinger et al 1997).
In this study, we compare the performance of the FENE-P model with FENE models with up to 20 beads. To overcome the computational expense of the multi-bead FENE model, we solve the equations along trajectories through a Newtonian DNS of isotropic turbulence. The resulting velocity information is stored and used to calculate the dynamics of all the polymer models. By comparing the performance of the various polymer models in the same flow field, we obtain a quantitative estimate of the error associated with the FENE-P model. Our approach is similar to the one taken by Zhou and Akhavan (2003) and others (Massah and Hanratty 1997, Massah et al 1993. However, those studies focused on inhomogeneous turbulent channel flow, where the statistics of the trajectory are a strong function of the initial condition. Obtaining statistical convergence for such a flow is difficult. In this study, we focus on statistically stationary isotropic turbulence, allowing us to perform spatial and temporal averages to obtain better statistical convergence. Additionally, we investigate the spatial resolution requirements for FENE-P simulations by considering the polymer stress variation along trajectories that have been pre-calculated to converge to a small volume in the flow domain. By examining the variation of the stress along neighboring points as a function of the size of the domain, we obtain guidance on the spatial resolution requirements for the FENE-P model as a function of the parameter values, varied over ranges of relevance to drag reduction.

Polymer models
Due to its complex molecular structure, a polymer molecule can have a broad range of length and time scales. Simulation techniques based on an atomistic description of the polymer, such as molecular dynamics (MD) and Monte-Carlo (MC) simulations, have been widely used to study polymers such as DNA and proteins. The atomistic description, however, is very expensive and is more suited to understanding the equilibrium configuration of the molecule and its chemical interactions with the solvent molecules, if present. MD and atomistic MC can span only picoto nano-seconds, and hence are not suitable for hydrodynamic studies.

FENE
To achieve the longer time scales required in hydrodynamic studies (i.e. micro-seconds to seconds), atomistic descriptions of the polymer molecule are typically 'coarse-grained' to consider only the longer relaxation times of the polymer. Excellent reviews of the literature on coarse-grained polymer models are provided by Bird et al (1987), Larson (1988) and Bird and Ottinger (1992). The models can be broadly categorized as bead-rod models such as 4 Kramers chain and bead-spring models such as the FENE model. These models use a bead to represent an inflexible segment of monomers (called a Kuhn length) in the polymer chain, and use rods or springs to model the interactions between the beads. With these simplifications, the degrees of freedom in the system are greatly reduced. The bombardment of the beads by the solvent molecules due to thermal fluctuations is modeled as a white noise stochastic process (Gardiner 1985, Ottinger 1996. Brownian dynamics (BD) simulations of this class of models are used to describe the extension and refolding of a polymer chain on hydrodynamic time scales. Examples of BD simulations of polymer models can be found in the work of Li et al (2000), Ghosh et al (2002b) and Terrapon et al (2004).
The governing equation for an N -bead FENE chain, denoted hereafter as FENE-N , is derived by applying a force balance on each of the beads (neglecting their accelerations). The force balance on the mth bead yields where F (m) is the force due to the spring connecting the mth and (m + 1)th beads, F (m) D is the force due to viscous drag and F (m) B is the Brownian force due to collisions of the solvent molecules. Note that in our convention, F (0) and F (N ) are zero, since the end beads have only a single spring attached to them. In the standard FENE formulation, the drag force is modeled as Stokes drag, the Brownian term is modeled by a Wiener process and the nonlinear spring force is modeled as a Warner (1972) spring. The resulting expressions are where u (m) is the undisturbed fluid velocity at the center of the mth bead, v (m) is the bead velocity, ζ ≡ 3π µ s d is the drag coefficient, µ s is the solvent viscosity, d is the effective bead diameter, k is the Boltzmann constant, T is the absolute temperature, W is a Wiener process, R (m) and R (m) ≡ |R (m) | are the separation vector and distance between the (m − 1)th and mth beads respectively and L max is the maximum extension of the spring. Notice that the restoring spring force diverges as R (m) → L max , preventing the separation distance from exceeding L max . Substituting these relationships into equation (1) and taking the difference of the equations for the (m + 1)th and mth beads yields The velocity difference between the two beads can be written as

5
If we invoke the 'locally linear flow assumption', the difference in fluid velocities is approximated by a Taylor series expansion to first order where ≡ ∇u is the velocity gradient tensor. Combining these relationships leads to the following stochastic differential equation for the separation vector It is common to introduce a spring relaxation time, λ(N ) ≡ ζ /4H , and equilibrium length scale, a ≡ √ kT /H , into the equation. If we define a normalized separation vector as r (m) ≡ R (m) /a, and a normalized maximum extension parameter as b(N ) ≡ H L 2 max /kT , then the above equation becomes The parameters λ(N ) and b(N ) are denoted as functions of the number of beads, N , because they are often adjusted to fix the corresponding properties of the chain, here denoted by λ and b, respectively. We discuss the relationships between the individual spring properties and the chain properties in section 2.3. The polymer stress exerted on the fluid is given by (Bird et al 1987) where · is the ensemble average over the configuration space of the polymer realizations, and n is the number density of the polymer molecules. A second quantity of interest is the elastic energy stored by the polymer. For the FENE model, this can be derived by integrating the Warner spring force law. The resulting relationship is where E is the stored elastic energy per unit volume, and P 0 is a reference energy. The above equations can be simulated using standard techniques for the white noise terms (Ottinger 1996).

FENE-P model
While the above FENE model provides the means to achieve time scales of relevance to hydrodynamics, it is still computationally expensive to integrate into a general fluid mechanics solver due to the need to establish a representative ensemble of FENE chains at each grid point. The CONNFFESSIT scheme is based on this approach (Laso and Ottinger 1993); however, calculations are generally limited to two-dimensional flows. A further simplification that is often applied in turbulence simulations is the closure of the FENE-2 (dumbbell) model introduced by Peterlin (1961). He recognized that if the restoring force were defined in terms of the average stretch (in place of the instantaneous value), the equations are closed. We begin with the 2-bead (dumbbell) form of equation (7) dr = r · dt + 1 λ dW − 1 2λ 2dW. Peterlin (1961) approximated this equation by where f (c) is the Peterlin function, defined as We define the conformation tensor as the mean of the outer product of the r vector with itself, i.e. C ≡ rr , and its closed transport equation can be derived from equation (10) by standard means (Pope 2000) dC where c = Tr(C). The expression for the polymer stress tensor now becomes where β ≡ µ S /µ 0 is the ratio of the solvent viscosity to the mixture viscosity in the limit of zero shear. β plays the role of a normalized polymer concentration. The polymer elastic energy for the FENE-P model is given by where again P 0 is the reference energy. Note that the normalization we have selected yields at equilibrium (i.e. in the absence of flow) Many researchers prefer to define the equilibrium condition as C eq = I, in which case the equations are modified (see appendix). 7

Parameters
The parameters that define the polymer chain are the ratio of the square of the maximum extension to its equilibrium extension, b, and the longest chain relaxation time, λ. For the FENE-N model, these molecular properties must be mapped into the properties for each individual spring, referred to respectively as b(N ) and λ(N ). Following convention, we express the relationship for the polymer relaxation time in terms of the coefficient K (N ) ≡ λ/λ(N ). Note that for a 2-bead chain (dumbbell) with a single spring, the chain properties are the same as the spring properties, i.e. b(2) = b and K (2) = 1. The Weissenberg number is a dimensionless ratio of the polymer chain relaxation time to a characteristic time of the flow. In isotropic turbulence, the Weissenberg number is generally defined as W e ≡ λ/τ η , where τ η ≡ (ν/ ) 1/2 is the small-scale (Kolmogorov) time scale. A third parameter is the polymer loading defined as n. In this study, we only consider one-way coupling, which is equivalent to taking the dilute limit, i.e. n → 0 or β → 1.
In order to make meaningful comparisons between the FENE-N and FENE-P models, we must define accurate relationships between the overall chain parameters and the individual spring parameters as a function of the number of beads, N . Indeed, we observe a strong sensitivity of the results to precisely how the parameter mapping is carried out.
We begin with the classical mappings used in viscoelastic laminar flows. According to the 'random walk theory' by de Gennes (1979), the equilibrium square end-to-end distance of the chain increases linearly with the number of springs. From the definition of b, we obtain This mapping is commonly used in the literature (Ghosh et al 2002a(Ghosh et al , 2002b. There exist several mappings for the spring relaxation parameter, K (N ). This is, in part, due to the fact that the ideal mapping is sensitive to the statistic chosen to be matched, as well as the nature of the flow. Rouse (1953) derived the following relationship from statistical thermodynamic considerations for weakly stretched chains (i.e. obeying a linear spring force law) Wiest and Tanner (1989) derive a separate expression based on a second-order Taylor series expansion for the zero-shear-rate material functions of the FENE-N model .
We shall refer to equations (16) and (17) as the 'Rouse' model and equations (16) and (18) as the Wiest and Tanner model, designated 'WT89'. The Rouse and WT89 mappings are optimized for polymers that are only moderately stretched. For highly stretched polymers, Wiest and Tanner (1989) showed that the elongational viscosity behaves asymptotically like 8 whereė is the elongational rate of the uni-axial extensional flow. Equation (19) implies that a separate constraint exists for highly stretched chains that is satisfied by neither the Rouse nor WT89 mappings. Zhou and Akhavan (2003) recognized this deficiency and suggested that K (N ) be defined by equation (17) and b(N ) be defined by forcing the right-hand side of equation (19) We shall refer to equations (17) and (20) as the Zhou and Akhavan or 'ZA03' mapping. Finally, we suggest a modification of ZA03. We assume b(N ) is given by the random walk theory shown in equation (16). Enforcing the right-hand side of equation (19) to be independent of N and K (2) = 1 implies Alternative approaches to derive equation (21) are shown by Jin (2007). We refer to the combination of equations (16) and (21) as the 'new' mapping. Below, we compare the performance of the four mappings: Rouse, WT89, ZA03 and new. As a preliminary test, we consider the performance of FENE-N and FENE-P in laminar flows. The first is uni-axial extension, defined by the uniform velocity gradient whereė is the extensional rate (strength) of the flow. Results for the FENE-N models are obtained by simulating an ensemble of 512 molecules and averaging. The analytical solution for the FENE-P model is where W e e ≡ėλ is the Weissenberg number for this flow, f is the Peterlin function (see equation (11)) and c ≡ Tr(C) is the real root of the following third-order polynomial circle. An ideal mapping in these coordinates is a horizontal line. We observe that the Rouse and WT89 mappings have significant variations in stress as a function of N , whereas ZA03 and the new mapping yield nearly identical straight lines. We attribute the improved performance of the latter two mappings to the fact that they both produce the same elongational viscosity for a highly stretched chain. With respect to end-to-end distance, the new mapping performs much better than ZA03, most likely due to the better performance of equation (16) for mapping b(N ).
Next, we consider laminar shear flow, defined by the velocity gradient where S is the shear rate. Once again, the FENE-N model is simulated numerically for an ensemble of 512 chains. The analytical solution for the FENE-P model for this flow is where W e S ≡ Sλ is the shear flow Weissenberg number and c ≡ Tr(C) is the real root of the third-order polynomial 2W e 2 S c 3 − 6bW e 2 S c 2 + (3b 2 + b 3 + 6b 2 W e 2 S )c − (2b 3 W e 2 S + 3b 3 ) = 0. Figure 2 shows the polymer stress and end-to-end distance for shear flow. The new mapping yields nearly a straight line (with the exception of a kink at N = 3 due to occasional folding of the chain). The next best performer is ZA03 followed by the Rouse and WT89 mappings. The end-to-end distance predictions of all the models are nearly identical. The error for FENE-P in this flow is significantly larger than for uni-axial extension. This is due to the complex dynamics observed with the FENE-N models. The FENE-N chains, stretched in the direction of the streamlines by the shear, occasionally tumble when the Brownian force causes the chain's orientation to pivot past the horizontal axis, so that the rear of the chain is exposed to the higherspeed fluid. During a tumble event, the chain completely retracts, before being re-stretched by the flow (click here to view an animation of a FENE-2 and FENE-20 molecule tumbling in a laminar shear flow). Tumbling in shear flow causes the FENE-N chain to sample a much broader region of the configurational phase space than the same chain will experience in uniaxial extensional flow. Consequently, the assumption made in deriving the FENE-P model is less valid, giving rise to the larger errors found in figure 2. Figures 1 and 2 show that the accuracy of the prediction of FENE-P depends upon the flow type. Flows that are more extensional in character are more accurately predicted than pure shear flow. The new and ZA03 mappings yield better overall performance than the Rouse and WT89 mappings, in the parameter range of interest to this study. However, the results thus far have been for steady state behavior in steady flows. Turbulent flows are inherently time dependent, and so the temporal response of the models to flow transients also will be important. Figure 3 demonstrates an important transient phenomenon known as the 'hysterisis effect'. The figure is obtained by turning on uni-axial extensional flow for a period of time, then turning it off and allowing the polymer to relax back to its equilibrium state. In the phase space showing the elongational viscosity, µ − 3µ s , versus the degree of extension, r 2 /b, we see the FENE-2 model exhibits a hysterisis that cannot be reproduced by the FENE-P model, highlighting a fundamental shortcoming of the FENE-P model in transient flows. Turbulence is inherently transient, hence it is necessary to test the accuracy of the predictions of the FENE-P model under turbulent flow conditions that are relevant for drag reduction investigations. Demonstration of the hysterisis effect. Transient normalized polymer elongational viscosity, µ − 3µ S versus the rate of extension r 2 /b. Uni-axial elongation is turned on first to stretch the polymer and then turned off to let it relax. The solid red line shows the phase space trajectory for a FENE dumbbell for W e e = 6, b = 50, and the black line with circles shows the same curve for the FENE-P model with the same parameters. Note that the FENE-P curve shows no evidence of the hysterisis found with the FENE-2 curve.

Comparison of polymer models in isotropic turbulence
Numerical simulations of FENE-2, FENE-5, FENE-10 and FENE-20 chains were performed along trajectories of forced Navier Stokes turbulence. The turbulence was generated by a standard Newtonian DNS code using 64 3 grid points and with forcing over the first two wavenumbers. The Reynolds number based on the Taylor microscale is R λ = 65 and κ max η = 1.2, where κ max is the maximum wavenumber of the simulation and η is the Kolmogorov length scale. Details of the numerical methods used in the simulations can be found elsewhere (Collins and Keswani 2004). The velocity gradient along 100 independent trajectories was spectrally interpolated, stored and later used to evolve the FENE-N and FENE-P models. Along each trajectory, an ensemble of 512 realizations of the model was calculated and averaged over. Figure 4 shows snapshots of the FENE-20 chain at different times along a single trajectory. Notice that at some times the molecule is stretched in a single direction (e.g. cases (a) and (b)), while at other times the molecule is more coiled (cases (c) and (d)). Coiling and uncoiling of the molecule occurs in response to fluctuations in the local strain rate and vorticity. Note that the caption to figure 4 contains a link to an animation of a FENE-2 and a FENE-20 chain along one of the isotropic turbulence trajectories. To evaluate the performance of the various FENE-N models and FENE-P, we must define the quantities to compare. Prediction of the polymer stress is critical, as this quantity enters directly into the momentum balance equation. The polymer stress tensor can be viewed in terms of its invariants and its orientation. To compare the prediction of the 'magnitude' of the polymer stress, we consider the first invariant or the trace of the polymer stress tensor defined as T [p] ≡ Tr(T [p] ). We consider FENE-20 to be the most accurate polymer simulation we performed, since it contains the largest number of degrees of freedom. Hence, we define the relative error for the FENE-N (N < 20) and FENE-P models as follows

12
where N (t) is the instantaneous percent error in the magnitude of the polymer stress. Figure 5 shows the time evolution of N (t) over a single trajectory for the four parameter mappings discussed earlier. The new and ZA03 mappings perform comparably and better than the Rouse and WT89 mappings at the relatively large values of Weissenberg number considered. We attribute this to the fact that the polymer remains highly stretched over most of the duration of the simulation. Thus mappings that preserve the extensional viscosity are more successful. We define the average error of the magnitude of the stress tensor, N , as the mean of | N (t)| over time and over the ensemble of 100 trajectories. We summarize the average errors in the trace of the polymer stress tensor in table 1. We see that the errors for FENE-P are on the order of 10% for the conditions considered. We note that instantaneous relative errors may be significantly larger, but this is less relevant for the application of drag reduction that is predominantly a global phenomenon.
To evaluate the errors in the 'orientation' of the polymer stress, we note that the polymer stress tensor is positive semi-definite, and hence its eigenvalues, λ [1] , λ [2] and λ [3] are real and positive, and the associated eigenvectors Λ [1] , Λ [2] and Λ [3] , form an orthogonal basis set. Note that these definitions apply equally to FENE-N chains, dumbbells and the FENE-P model. Without loss of generality, we adopt the convention λ [1] > λ [2] > λ [3] . Orientation can be measured by comparing the directions of the eigenvectors. We are mainly concerned with Λ [1] , since this provides the orientation of the dominant force exerted on the fluid by the polymer.      We define an instantaneous angle between various FENE-N chains and FENE-20 as (23) Figure 6 shows a comparison of θ N (t) in radians for the four mappings. Here, we see that all of the mappings yield nearly identical, very accurate predictions of the stress orientation. As the polymers are highly stretched, and only rarely recoil, the orientation is well captured by the FENE-P model. This is evident in the animation linked to the caption of figure 4 showing that the FENE-2 and FENE-20 models are in good agreement with each other over most of the simulation. The largest errors in orientation occurring at t/τ η ∼ 60 and 80 coincide with the largest errors in the stress magnitude, N (t), highlighting the sensitive coupling between orientation and stress magnitude. We define the average angle, θ N , as the average of θ N (t) over time and the 100 trajectories, and report these values in table 2. You can see that the average orientation errors are extremely small, even for the FENE-P model. We can gain some insight into why the angles of the polymer stress tensor are so well predicted by the FENE-P model by looking more closely at how the polymer is aligning in the flow. Figure 7 shows the probability density function of the cosine of angles between vorticity (left) and the principal eigenvector of the polymer stress tensor (right) with the eigenvectors of the rate-of-strain tensor (ordered by eigenvalues, largest to smallest), and with a material line obeying dr = r · dt.
The plot on the left shows the tendency of vorticity to align with the second eigenvector of the rate-of-strain tensor, which has been discussed earlier (Ashurst et al 1987). There is a somewhat similar tendency for the polymer stress tensor eigenvalue, although the alignment with both the first and second rate-of-strain eigenvectors shows a positive correlation. There is a clear negative correlation with the third rate-of-strain eigenvector (the compressional axis). What is most interesting is how closely aligned the polymer stress tensor eigenvector is with the material line. If you neglect the Brownian (equilibrium) terms in the FENE-N (FENE-P) equation, which is a valid assumption at high Weissenberg numbers, the equation for the orientation of the polymer reduces to a form that is nearly the same as the equation for a material line. The only exception is due to the restoring force that acts along the polymer axis for a highly stretched chain, and hence does not affect the orientation of the polymer.
We conclude that FENE-P provides a very good description of the polymer stress tensor, both in terms of its magnitude and orientation. This conclusion is more optimistic than the one drawn by Zhou and Akhavan (2003), who focused more on the instantaneous error than the average error. We likewise observe instantaneous errors in the polymer stress magnitude in excess of 100%, but for a comparably small percentage of the time. We believe the average error is more indicative of the overall performance of the model in terms of its prediction of average properties of the turbulence-polymer interaction such as percent drag reduction.

Linearity and resolution
In this section, we investigate the spatial resolution requirement for the FENE-P model using Lagrangian simulations along trajectories that have been pre-calculated to terminate on subgrid mesh points surrounding one of the grid points (see figure 8). The study is based on a 64 3 DNS of stationary isotropic Newtonian turbulence (R λ = 55, κ max η = 1.5). Identical Lagrangian simulations are performed on 64 evenly spaced grid points throughout the simulation domain. For each grid point, we identify 7 3 = 343 subgrid points that uniformly surround a given grid point. We define the ratio of the grid separation distance to the subgrid separation distance as ς ≡ x/ x and set ς = 1, 20, 100, 10 4 and 10 6 in this study, where ς = 1 corresponds to the original DNS grid resolution. It would be difficult to identify the starting location for all 343 trajectories that fall precisely at the nodes of the sub-grid on the final time step. Instead, we take advantage of the stationary nature of the turbulence and initialize the trajectories at the nodal points and simulate them forward in time, storing the velocity gradient along each of the trajectories. We then invert the trajectories so as to attain velocity gradient statistics along Figure 9. A demonstration of trajectories coming together on to the sub-grid surrounding a grid point. The trajectories begin from seemingly random spatial locations, but as time goes on they come together and finally land precisely on the sub-grid points surrounding one grid point. Click here to view an animation of the converging trajectories.
trajectories that end at the sub-grid nodes. Note that this is not the same as truly inverting time, since the velocities (and velocity gradients) would change sign (inverting extensional and compressional axes for strain, etc). Here, we simply have re-ordered the velocity gradients with the correct signs based on stationarity. We believe the inverted trajectories provide at least qualitative information on the sub-grid convergence. Figure 9 shows an example of the trajectories used.
A continuous field is considered well resolved by a grid when the field values between adjacent points vary nearly linearly. Conversely, if a linear interpolation is unable to approximate the variables between grid points, the grid spacing is insufficient to resolve the field. A tensor field is 'linear' if it can be well represented by Given at least 27 grid values of i j (x), we can perform a least squares fit and obtain the coefficient of determination, 2 ∈ [0, 1], as a measure of the degree of linearity (a linear fit corresponds to 2 ≈ 1). The mathematical definition of the coefficient of determination can be found in the book by Abraham and Ledolter (2005). The linearity measure for the velocity gradient itself is shown in figure 10. As the velocity gradient is well resolved (k max η = 1.5), it is not surprising that we find that 2 approaches unity for all values of ς . We observe the period of time over which 2 ≈ 1 increases with increasing ς . This is due to the fact that the period of time during which the converging trajectories are in close proximity to each other increases with increasing ς .
Results for the linearity of the polymer stress computed from the FENE-P simulations along the trajectories are shown in figure 11. We notice that the maximum value of 2 (corresponding to t/τ η = 0) does not approach unity for the smaller values of ς , suggesting that the original grid spacing (corresponding to ς = 1) cannot fully resolve the polymer stress field at the moderate-to-high Weissenberg numbers considered in this study. The cases corresponding to ς = 20 and 100 achieve reasonably high values of 2 , while the finer-grained cases yield values very close to unity. As noted above, as ς increases, the period of time during which the converging trajectories are in close proximity to each other increases, enabling the initially uncorrelated values of polymer stress to correlate, causing 2 to increase. The final value of 2 , therefore, is related to the correlation time of the trajectories as they approach each other.
Another observation from figure 11 is that when the Weissenberg numbers are large, the linearity is not very sensitive to its value. This is encouraging because it implies that the required spatial resolution for a DNS eventually becomes independent of the Weissenberg number.
These results suggest that a much finer grid is required to fully resolve the polymer field than is required to resolve the velocity and pressure fields. This is qualitatively similar to the problem of shock waves in gas dynamic flows, and results from the hyperbolic character of the FENE-P model. Full resolution of a gas dynamic shock is not practical, and hence numerical schemes such as Lax-Friedrichs, Lax Wendroff and the Godunov methods have proven successful at describing the shock wave (e.g. its strength, propagation speed, etc) without fully resolving it. Motivated by these results, Vaithianathan et al (2006) adapted the secondorder hyperbolic solver of Kurganov and Tadmor (2000) to advance the conformation tensor equations. The advantage of this flux limiter-based scheme is that it adjusts in the vicinity of shocks so that the bounds on the eigenvalues cannot be violated, eliminating the instabilities that can arise in these types of calculations (Vaithianathan and Collins 2003), without introducing a global stress diffusivity (Sureshkumar and Beris 1995).

Conclusions
In this study, we examined the accuracy of the FENE-P model for polymers in turbulence by computing the polymer stress tensor along individual trajectories in an isotropic DNS of Newtonian turbulence. Parameters were set to values of relevance to polymer drag reduction, namely we focused on relatively high values of the Weissenberg number, since this represents the class of polymers most effective for this application. We benchmarked FENE-P by comparing its performance to FENE-20. We also investigated the performance of FENE-10, FENE-5 and FENE-2. We find errors in the magnitude of the polymer stress of approximately 10-15% over all conditions studied. We also investigated the orientation statistics of the polymer stress by comparing the direction of the dominant eigenvector (eigenvector associated with the largest eigenvalue) of the polymer stress tensor. We find this angle to be very small (less than one degree on average). This is due to the fact that the polymer chains are highly stretched and hence behave nearly like a material line in the flow, independent of the number of beads in the chain.
The second aspect of the model we investigated is the requirement for spatial resolution of the FENE-P model relative to the requirement for the flow. We note that the equation for the conformation tensor C is hyperbolic, and hence capable of producing discontinuities in the magnitude and orientation of the tensor. To investigate this, we performed Lagrangian simulations of FENE-P along trajectories that were pre-calculated to converge to a subgrid of points surrounding a particular node of the computation mesh. These trajectories were generated by solving the inverse problem, namely the divergence of trajectories initially on the subgrid, 20 and then inverting the temporal order of the trajectories. The results indicate that the polymer stress tensor can have much sharper gradients than the underlying flow, but not necessarily discontinuities. We trace the smoothing to the relaxation term in equation (12), which causes the polymer stress in neighboring points to correlate, even though their initial flow histories were uncorrelated. Nevertheless, the sharp features found in the polymer stress suggest that it cannot be easily resolved in a three-dimensional DNS, motivating the use of hyperbolic solvers for these equations (Vaithianathan et al 2006).
Readers interested in performing numerical simulations of the FENE-N or FENE-P models can download our open source code called QPolymer (http://sourceforge.net/projects/ qpolymer/). The code allows you to simulate an arbitrary FENE-N chain along a model turbulent flow trajectory. Details of how to install and use the code can be found on the web site.
The stress tensor is given by (A.6) and the polymer elastic energy by where P 0 is the reference energy. The correct formulation is either the equations in the main text or the equations in the appendix. We note there are several formulations in the literature that improperly combine these two normalizations (e.g. using the Peterlin function from one formulation and the governing equations from the other). The errors are small for b 3, but can be significant for smaller values of b.