Elastic waves and transition to elastic turbulence in a two-dimensional viscoelastic Kolmogorov flow
We investigate the dynamics of the two-dimensional periodic Kolmogorov flow of a viscoelastic fluid, described by the Oldroyd-B model, by means of direct numerical simulations. Above a critical Weissenberg number the flow displays a transition from stationary to randomly fluctuating states, via periodic ones. The increasing complexity of the flow in both time and space at progressively higher values of elasticity accompanies the establishment of mixing features. The peculiar dynamical behavior observed in the simulations is found to be related to the appearance of filamental propagating patterns, which develop even in the limit of very small inertial non-linearities, thanks to the feedback of elastic forces on the flow.
It is well known that the addition of small amounts of long chain polymers can have a strong impact on flowing fluids’ properties. One of the most remarkable effects of highly viscous polymer solutions which has been recently observed in experiments is the development of an “elastic turbulence” regime in the limit of small Reynolds number and strong elasticity (1). The flow in this regime displays a rich dynamics, of the type usually encountered in high Reynolds number turbulent flow in Newtonian fluids. Typical features of turbulent flows (e.g., broad range of active scales, growth of flow resistance) are observed in the experiments by Steinberg and collaborators (1); (2), even at low velocity and high viscosity, i.e. in the limit of vanishing Reynolds number. As a consequence of turbulent-like motion, elastic turbulence has been proposed as an efficient technique for mixing in very low Reynolds flows, such as in microchannel flows (3); (4); (5). Despite its great technological interest, elastic turbulence is still only partially understood from a theoretical point of view. Presently available theoretical predictions are based on simplified versions of viscoelastic models and on the analogy with magneto-hydrodynamics (MHD) equations (6); (7). In a previous work (8) we showed, by means of numerical simulations, that the basic phenomenology of elastic turbulence is reproduced by a simple model of polymeric fluid and in an idealized geometry, namely a two-dimensional periodic shear flow without boundaries. Notwithstanding quantitative differences, most observed features have a strong qualitative resemblance with experimental results.
The mechanism through which polymer molecules can influence properties of a flow is their extreme extensibility. Polymers, typically composed by a large number of monomers, at equilibrium are coiled in a ball of radius . In presence of a non-homogeneous flow, the molecule is deformed in an elongated structure characterized by its end-to-end distance which can be much larger than . The deformation of molecules is the result of the competition between the stretching induced by velocity gradients and the entropic relaxation of polymers to their equilibrium configuration. Experiments with DNA molecules (9); (10) show that this relaxation is linear, provided that the elongation is small compared with the maximal extension , and therefore can be characterized by a typical relaxation time (11). Besides the Reynolds number , commonly used in hydrodynamics, in the study of flowing polymer solutions, it is useful to introduce the Weissenberg number . This is defined as the product of and the characteristic deformation rate, and it measures the relative importance between polymer relaxation and stretching. When relaxation is fast compared to the stretching time and polymers remain in the coiled state. For , on the contrary, polymers are substantially elongated. The transition point is called the coil-stretch transition and occurs at . For Weissenberg numbers much larger than unity, polymers can provide a significant feedback effect on the dynamical behavior of the fluid they are suspended into, so that they behave as active objects capable of modifying the properties of the flow.
Because of the feedback effect, above the coil-stretch transition, a laminar flow can become unstable even in the small limit of vanishing inertial non-linearities. Purely elastic instabilities have been studied in different flow configurations. In a curvilinear shear flow as, e.g., Couette flow between rotating cylinders, the driving force of the elastic instability is the “hoop stress” (see, for instance, (12) page 62), a volume force in the direction of curvature due to the first normal stress difference. The mechanism for this instability was first proposed in (13) and experimentally verified in (14). In the framework of flows with curvilinear streamlines, a dimensionless criterion for the critical onset conditions, which applies to various flow geometries, was established in (15); (16). The role of flow-microstructure coupling in pattern formation in the Taylor-Couette configuration was also examined in numerical simulations of the FENE-P (finitely extensible non-linear elastic-Peterlin) model (17). On the other hand, purely elastic instabilities in a parallel flow with rectilinear streamlines were obtained for the first time in a Kolmogorov flow of a viscoelastic fluid described by the linear Oldroyd-B model (18). Also for a more realistic shear flow of a viscoelastic fluid in a straight channel, there exist theoretical results exhibiting the elastic instability (19), although they are still not experimentally verified. Elastically induced pattern formation and production of fluid mixing were also studied in extensional flows, in experiments (20), and numerical simulations of a Stokesian linear viscoelastic model (21).
In this paper, we study the dynamics of dilute polymer solutions by means of direct numerical simulations of a general linear viscoelastic model, in a flow configuration corresponding to the two-dimensional Kolmogorov shear flow. Our simulations for values of the stability parameters ( and ) close to the elastic instability threshold, exhibit a transition to new intriguing flow states, in the form of “elastic waves”, driving the establishment of mixing features. In particular, at growing the elasticity of the solution, after a first transition to steady uniformly traveling patterns with different symmetry with respect to the reference laminar fixed point, we detect the appearance of temporal oscillations and, eventually, the onset of a strongly non-linear multiscale state disordered in both space and time.
The article is organized as follows: in the next section we introduce the viscoelastic model. In section III.1 we present numerical results concerning the transition to mixing states in terms of an analysis of the behavior of global observables. In section III.2 we show the results of numerical measurements concerning the spatiotemporal structure of the patterns emerging above the onset of the instability. Finally, section IV summarizes our conclusions.
To describe the dynamics of a dilute polymer solution we adopt the well known linear Oldroyd-B model (12)
where is the incompressible velocity field, the symmetric positive definite matrix representing the normalized conformation tensor of polymer molecules and is the unit tensor. The trace is a measure of local polymer elongation. The solvent viscosity is denoted by and is the zero-shear contribution of polymers to the total solution viscosity and is proportional to the polymer concentration. In absence of flow, , polymers relax to the equilibrium configuration and . The extra stress term in (1) takes into account the elastic forces of polymers and provides a feedback mechanism on the flow.
In order to concentrate on the intrinsic dynamics of the viscoelastic solution we use the simple geometrical setup of the periodic Kolmogorov flow in two dimensions (22). With the forcing , the system of equations (1-2) has a laminar Kolmogorov fixed point given by
with (18). The laminar flow fixes a characteristic scale , velocity and time . In terms of these variables, one can define the Reynolds number as and the Weissenberg number as . The ratio of these numbers defines the elasticity of the flow . Notice these are a-priori estimates of the non-dimensional parameters; in the following we will refer to them as and . In the presence of velocity fluctuations superimposed to a non-zero mean flow , it is also possible to compute a-posteriori non-dimensional parameters as and , once the average velocity has been computed; of course, in general .
It is well known that the Kolmogorov flow displays instability with respect to large-scale perturbations, i.e. with wavelength much larger than . In the Newtonian case, the instability arises at (23). The presence of polymers changes the stability diagram of laminar flows (24); (25) and can induce elastic instabilities which are not present in the Newtonian limit (13); (26); (18); (14), even in the case of periodic flows (27). In this respect, the Kolmogorov flow is no exception, and recent analytical and numerical investigations have offered the complete instability diagram in the (, ) plane (18). For the purpose of the present work, we just have to recall that linear stability analysis shows that for sufficiently large values of elasticity, the Kolmogorov flow displays purely elastic instabilities, even at vanishing Reynolds number (see Fig. 1 of (18)). Note that in this case the direction of the most unstable mode is at a small angle with respect to the basic flow, at variance with the purely hydrodynamic transverse instabilities. In the Newtonian case, non-transverse instabilities are observed for less symmetric basic flows only (28). We also remark that the fact that the basic flow has rectilinear streamlines does not exclude the onset of the elastic instability (19). Above the elastic instability the flow can develop a disordered secondary flow which persists in the limit of vanishing Reynolds number and eventually leads to the elastic turbulence regime (29).
Iii Analysis and results
The equations of motion (1)-(2) are integrated by means of a pseudo-spectral method implemented on a two-dimensional grid of size with periodic boundary conditions at resolution . It is well known that the integration of viscoelastic models is limited by numerical instabilities associated with the loss of positiveness of the conformation tensor (30). These instabilities are particularly important at high elasticity and limit the possibility to investigate the elastic turbulent regime by direct implementation of equations (1)-(2). To overcome this problem, we have implemented an algorithm based on a Cholesky decomposition of the conformation matrix that ensures symmetry and positive definiteness (31). While this allows us to safely reach high elasticity it has a price to pay in terms of limited resolution and large computing time. In the following, we will present the numerical results obtained using this approach.
iii.1 Transition to mixing states
In order to investigate the transition to mixing states, we chose to work in a region of parameters close to the onset of the purely elastic instability. Here we fix and increase the Weissenberg number from to . The initial condition is obtained by adding a small random perturbation to the fixed point solution (3) and the system is evolved in time until a (statistically) steady state is reached.
We start the analysis by considering the behavior of global quantities, such as the kinetic energy and the average square polymer elongation . The total energy of the system is the sum of kinetic and elastic contributions, and its balance equation, from (1)-(2) is:
where is the energy input, is the viscous dissipation, and the last term represents elastic dissipation. The balance of forcing and dissipation produces a steady state, in which the total energy is constant, at least in a statistical sense (i.e. on time scales larger than the typical duration of fluctuations).
The mean values of and as a function of the Weissenberg number are presented in Fig. 1. Here and in the following the Weissenberg number and the Reynolds number are defined in terms of the amplitude of the mean velocity, which is measured from velocity profiles. Indeed, a remarkable feature of the Kolmogorov flow is that even in the turbulent regime the mean velocity and conformation tensor are accurately described by sinusoidal profiles (32), though with different amplitudes with respect to the laminar fixed point.
As it is shown in Fig. 1, the system does not depart from the laminar fixed point up to which is approximately the critical obtained in (18). For higher values of the Weissenberg number, a depletion in the kinetic energy is observed, accompanied by a growth in the trace of the conformation tensor, larger than its laminar prediction . Therefore, polymers drain energy from the mean flow in order to elongate and, once sufficiently stretched they appreciably start to back-react on the flow itself, through the coupling term in (2). The increase with respect to the laminar growth indicates this coupling is very efficient in sustaining the stretching of polymers.
We observe here that, due to the present flow geometry, one of the diagonal elements of the conformation tensor is much larger than the other, namely . This implies that the trace of the tensor is almost entirely given by , and the same is true for the first normal stress difference , so that . Together with the growth of the trace of the conformation tensor, the growth of the first normal stress difference is a clear signature of the elasticity of the flow. An inspection at the behavior of the components of the conformation tensor, Fig. 2, indeed reveals that and, in turn, . In particular we find independently of , and for . Although and always remain small compared to , their behavior with Weissenberg number clearly shows a transition around , indicating the establishment of a nonzero transversal force, growing with , potentially capable of sustaining a secondary flow with nonzero component in the direction perpendicular to that of the base flow.
Let us now discuss the behavior of temporal fluctuations of the global quantities (kinetic energy) and (trace of the conformation tensor), when the Weissenberg number is increased, which provides further information on the change of dynamical regime. In Fig. 3 we report and versus time, for several values of , above the elastic instability threshold . The picture emerging from this figure is that of a transition from constant to randomly fluctuating states, via periodic ones.
Below , and are constant, although at different values with respect to the laminar fixed point. We remark that the constant behavior does not necessarily mean the flow is time independent, as it is also compatible with the presence of steady, uniformly traveling patterns of the form propagating with velocity .
When the kinetic energy and the squared polymer elongation oscillate around their average values. To assess the frequency content of the above fluctuating temporal signals, we performed a Fourier analysis, computing the power spectrum
where is an interval of time in the statistically steady state; the results are shown in Fig. 4 (similar results are obtained for ). This procedure allowed to identify two transitions: at an explicit time dependency appears, signaled by a single frequency peak in the spectra of and . In this range of parameters the system is in a periodic state, characterized by global oscillations in time. The peak frequency is found to be slightly larger than the frequency associated to the advective time scale, , and smaller than that associated to the polymer relaxation time (which is in turn much smaller than the frequency related to the inverse velocity gradient of the mean flow): . The time dependency of the kinetic energy implies that of the velocity field . It is known from the theory of dynamical systems that this is a sufficient condition for the appearance of Lagrangian chaos in a incompressible flow and the chaoticity of Lagrangian trajectories implies the establishment of mixing features (33). In previous numerical simulations (8), we measured the Lagrangian Lyapunov exponent , defined as the mean rate of separation of two infinitesimally close particles transported by the flow. This quantity, which gives an estimation of the inverse mixing time, was indeed found to be positive in a corresponding range of values, implying the onset of chaotic (Lagrangian) dynamics.
The present picture persists at increasing Weissenberg numbers, until at a second frequency, associated to a slower time scale, shows up. For several frequencies are present, indicating more disordered oscillations. In this regime, the spectral amplitudes seem to approximately follow an exponential decay.
In this latter case, the random fluctuations in time of the global quantities point at the onset of irregular dynamics associated to a chaotic, mixing, multiscale state. A fixed time snapshot of the vorticity field at is shown in Fig. 5 where an irregular, disordered pattern is clearly observable superimposed to the basic periodic flow.
The statistical characterization of this turbulent-like state is given by the spectrum of velocity fluctuations in the wavenumber domain. As shown in Fig. 6, the spectrum displays a power-law behavior , with exponent , reflecting the activity of a whole interval of spatial modes. The fast decay of the spectrum, tells us that the velocity and velocity gradients are primarily determined by the integral scale, i.e. on the system size. This implies that elastic turbulence corresponds to a temporally random, spatially smooth flow, dominated by strongly non-linear interactions between few large scale modes.
iii.2 Wavy patterns
As discussed in the previous section, for values of the Weissenberg number larger than , we observe a transition to non stationary states. We now focus on features of such states in space and in time.
For a typical filamental pattern emerges in both velocity (or vorticity) and fields. Snapshots of these fields in the plane are shown in Fig. 7. These filamental structures are found to perform a wavy motion along the horizontal direction of the mean flow. The possibility of observing “elastic waves” in polymer solutions was theoretically predicted within a simplified uniaxial elastic model (7) which has strong formal analogies with MHD equations, but their experimental observation is still lacking.
In a narrow range of the patterns are found to simply translate without changing their form. Already at temporal oscillations of pattern shapes are manifest, particularly on filament tails and in regions between them. Nevertheless, wavy motion is still clearly detectable. At larger Weissenberg number, filamental patterns are significantly less stable. Indeed they are observed to form and disrupt in the course of time, as they travel in both positive and negative horizontal directions.
Since patterns described so far travel along the longitudinal direction, it is natural to wonder what their typical speed is, how this possibly varies with the transversal coordinate, and how these quantities change when the elasticity of the solution is increased. We measured the speed of these elastic waves by tracking peaks of intensity in the fields , (), and following their motion along the direction (Fig. 8).
The results relative to the cases where coherent motion of structures was detected are reported in Table 1 and Fig. 9. In all the considered cases we found an average wave speed intensity close to half the peak fluid velocity: . Concerning spatial variations of , these show little dispersion around the mean at , as it can be seen also in the plot (Fig. 9) of the profiles along the transversal coordinate , which are essentially constant. The speed of different parts of a wave, for instance its tip and its tail, tracked by cross-sections at different coordinates, move at approximately the same speed, pointing to rigid motion of the structure, which travels substantially undeformed. We remark these dynamics are compatible with the constancy of the global quantities shown in Fig. 3, due to the fact that and are integrals over the space domain by definition, and to the presence of periodic boundary conditions. When the Weissenberg number is increased ( in Fig. 9a), wave speed profiles become less constant, indicating more complex motion. The larger scatter of values reflects the larger uncertainty in the measure of the propagation speed of patterns, which are now found to swing at -locations corresponding to their tails, and the presence of secondary less steady filamental structures in regions between the better defined primary waves. At even larger values of the Weissenberg number (), waves interact considerably more, starting to move in both directions and to break and reform after they “collide”; the distribution of wave speeds now exhibits significantly more variability along the direction. In fact, at the present Reynolds number (), this is the largest at which we succeeded in tracking waves. Above this value, patterns form more easily in various parts of the spatial domain, but interactions among them dominate the dynamics, producing a disordered flow in both space and time.
At the smallest value of we find that the dependence of the mean velocity on is very weak (see Fig. 9b), confirming the elastic nature of these instabilities.
In order to get some insight into the physical mechanism producing filamental patterns, we now consider the spatial structure of the components , of the velocity field, for and . The velocity field can be thought as the sum of the base flow and a secondary flow: , ; at the laminar fixed point, of course , . Here we are especially concerned with the establishment of a nonzero transversal velocity component and the breaking of translational invariance along the longitudinal direction, above the elastic instability threshold. Snapshots at fixed time of the fields and , computed from the vorticity field shown in Fig. 7a, are presented in Fig. 10a, 10b.
In the horizontal component of velocity , the periodicity of the base flow is clearly visible. The transversal component , on the other hand, displays localized quadrupolar patterns of positive and negative values, on a background of zero velocity. Just to fix the ideas, consider, for instance, the pattern centered around , being the height of the computational domain, on the right part of Fig. 10a. Notice that the position of this pattern matches exactly that of the tip of the most clearly visible wave in the vorticity field (Fig. 7a). Indeed, this structure is found to be stable and to travel along the horizontal direction, following the propagation of the vorticity filament. The lobes of these quadrupolar patterns (with white, black indicating positive, negative values respectively), appear in correspondence to locations of strong primary strain, where gradients of the base flow attain maximal values in magnitude. In these regions polymers are significantly stretched, by the terms and in eq. (2). In Fig. 10c a snapshot, at the same time, of the trace of the conformation tensor, which measures the square elongation of polymers, is reported. Strong elongations appear in thin filamental structures localized around coordinates where changes sign or, equivalently, where is maximum. As a consequence, it is in this regions that they can exert a relevant feedback on the flow, through the elastic stress term .
In other words, the transverse flow develops where the feedback of polymers on the velocity dynamics is strongest. Finally, the establishment of the transversal component of the flow is accompanied by a reduction of the longitudinal one, due to incompressibility (Fig. 10b). The result of this whole process is the formation of recirculation zones, which are evident in a plot of the corresponding stream function (such that , ) (Fig. 11). These in turn result to be organized in an array of counter-rotating vortices of equal strength, with a basic periodicity of in the longitudinal direction, which propagates horizontally.
Concerning the feedback of polymers on the flow field, the transversal velocity is found to be sustained by a persistent correlation between the transversal components of the velocity field and the elastic body force , confirming what anticipated in the preceding discussion. The profiles of these quantities are shown in Fig. 12, at , that is immediately on the left of the center of the quadrupolar pattern discussed above. Though this figure refers to a single realization at fixed time, this positive correlation is a feature persisting over time. In Fig. 12 it is also shown the profile of , which is the sign of the elastic force’s power. Though this fluctuates wildly, it is clearly positive in finite size regions located around , corresponding to the extension of the quadrupolar pattern on the right in fig 10a. This observation further supports the role of elastic contributions in the establishment of the transversal flow, and hence in the appearance of filamental structures.
To conclude, the picture emerging from this description is that filamental patterns are associated to the formation of arrays of counter-rotating vortices, sustained by the dynamical coupling between velocity and localized structures where polymers are strongly elongated. Indeed this is probably the simplest way to break the translational symmetry of the velocity field along the horizontal direction. The heuristic argument developed here is summarized in the scheme in Fig. 13. As a final remark, it may be interesting to observe that this allows to understand the direction of wave propagation. This results from the direction of the base velocity field, from which vortices of intensity of order form. Due to the periodicity of the base flow, both positive and negative horizontal directions of propagation are possible. The selected direction then only depends on where the positive correlation between a random perturbation of and the elastic body force happens to be established. As shown in Fig. 14, at increasing the elasticity of the solution, the probability to find positive correlations gets larger, which implies the possibility of forming filaments in a larger fraction of the domain. This reflects the experimental observation that when is sufficiently large, filaments’ propagation may happen in both directions.
Iv Summary and conclusions
In this work we have studied the dynamics of a two-dimensional parallel flow, namely the Kolmogorov shear flow, of a viscoelastic fluid, close to the elastic instability threshold predicted in (18) within the framework of linear stability analysis. By means of direct numerical simulations of the Oldroyd-B model, we investigated the destabilization of this flow induced by the elastic forces associated to the dynamics of polymer molecules in the solution. Above a critical Weissenberg number a transition to new dynamical states was observed. In this regime of large elasticity, the system was found to be characterized by the presence of filamental structures propagating in the direction of the mean flow and global oscillations introducing temporal dependencies in the observables. This phenomenon allows the establishment of mixing features. Close to the critical Weissenberg number the structures can be described as traveling waves associated to steady patterns of counter-rotating vortices sustained by the coupling between elastic forces and velocity field. The onset of wavy motion of filaments was found to occur even in the limit of very small Reynolds numbers, typically for , at .
Increasing the Weissenberg number, the system evolves towards more complex dynamics, where filaments are observed to continuously form and break as a consequence of interactions among them. As a consequence, several time scales and finer spatial scales result to be involved in the dynamical behavior, which displays strongly non-linear features resembling those of a turbulent flow. For instance, when (and ) the power spectrum of velocity fluctuations is given by a power-law of exponent larger than , in close agreement with experimental findings (1).
The present results refer to a simplified viscoelastic model in an idealized geometrical configuration. Nevertheless, this simple setup proved useful to detect the appearance of coherent structures in the form of “elastic waves”, driving the transition towards mixing states, and to reproduce the basic phenomenology observed in experiments. We hope these findings may help in the understanding of the phenomenon of elastic turbulence in more realistic configurations and in the comprehension of the elastic instability mechanism in parallel flows.
We are grateful to S. Rafaï for interesting discussions and a critical reading of the manuscript. S.B. acknowledges financial support from CNRS.
- A. Groisman and V. Steinberg, Nature 405, 53 (2000).
- T. Burghelea, E. Segre and V. Steinberg, Phys. Fluids 19, 053104 (2007).
- A. Groisman and V. Steinberg, Nature 410, 905 (2001).
- T. Burghelea, E. Segre, I. Bar-Joseph, A. Groisman and V. Steinberg, Phys. Rev. E 69, 066305 (2004).
- P.E. Arratia, G.A. Voth and J.P. Gollub, Phys. Fluids 17 053102 (2005).
- E. Balkovsky, A. Fouxon and V. Lebedev, Phys. Rev. E 64, 056301 (2001).
- A. Fouxon and V. Lebedev, Phys. Fluids 15, 2060 (2003).
- S. Berti, A. Bistagnino, G. Boffetta, A. Celani, S. Musacchio, Phys. Rev. E 77, 055306(R) (2008).
- T.T. Perkins, D.E. Smith and S.Chu, Science 264, 822 (1994).
- S.R. Quake S.R., H. Babcock and S. Chu, Nature 388, 151 (1997).
- E.J. Hinch, Phys. Fluids 20, S22 (1997).
- R.B. Bird, C.F. Curtiss, R.C. Armstrong, and O. Hassager, Dynamics of polymeric fluids, Wiley, New York (1987).
- R.G. Larson, E.S.G. Shaqfeh and S.J. Muller, J. Fluid Mech. 218, 573 (1990).
- A. Groisman and V. Steinberg, Phys. Fluids 10, 2451 (1998).
- P. Pakdel and G.H. McKinley, Phys. Rev. Lett. 77, 2459 (1996).
- G.H. McKinley, P. Pakdel, A. Öztekin, J. Non-Newtonian Fluid Mech. 67, 19 (1996).
- D.G. Thomas, R. Sureshkumar and B. Khomami, Phys. Rev. Lett. 97 054501 (2006).
- G. Boffetta, A. Celani, A. Mazzino, A. Puliafito and M. Vergassola, J. Fluid Mech. 523, 161 (2005).
- A.N. Morozov and W. van Saarloos, Phys. Reports 447, 112 (2007).
- P.E. Arratia, C.C. Thomas, J. Diorio and J.P. Gollub, Phys. Rev. Lett. 96, 144502 (2006).
- B. Thomases, M. Shelley, Phys. Rev. Lett. 103, 094501 (2009)
- V.I. Arnold and L. Meshalkin, Uspekhi Mat. Nauk 15, 247 (1960).
- L. Meshalkin and Y.G. Sinai, J. Appl. Math. Mech. 25, 1700 (1961).
- R.G. Larson, Rheol. Acta 31, 213 (1992).
- A. Groisman and V. Steinberg, Phys. Rev. Lett. 77, 1480 (1996).
- E.S.G. Shaqfeh, Annu. Rev. Fluid Mech. 28, 129 (1999).
- K. Arora, R. Sureshkumar and B. Khomami, J. Non-Newtonian Fluid Mech. 108, 209 (2002).
- B. Dubrulle and U. Frisch, Phys. Rev. A 43, 5355 (1991).
- A. Groisman and V. Steinberg, New J. Phys. 6, 29 (2004).
- R. Sureshkumar and A.N. Beris, J. Non-Newtonian Fluid Mech. 60, 53 (1995).
- T. Vaithianathan and L.R. Collins, J. Comp. Phys. 187, 1 (2003).
- G. Boffetta, A. Celani and A. Mazzino, Phys. Rev. E 71, 036307 (2005).
- M. Cencini, F. Cecconi and A. Vulpiani, Chaos. From simple models to complex systems, World Scientific, Singapore (2010).