Turbulence on a Desktop

This column is the last one that we will do as editors. We have been editing the Computer Simulations column continuously since our first column in the July/August 1989 issue of Computers in Physics. We wish to thank all of the authors who have contributed to the column and patiently suffered through our suggestions and revisions, and the editors of Computers in Physics and Computing in Science & Engineering for their help and support. We will really miss editing the column, in part because it gave us a vehicle for learning new physics and algorithms, and in part because we could ask the authors clarifying questions that would have been too embarrassing to ask in other contexts. Our reason for leaving is that as of 1 July 2001, Jan Tobochnik will become the new editor of the American Journal of Physics and Harvey Gould will become the new associate editor. We look forward to the challenges of editing the American Journal of Physics, and, as usual, we look forward to your contributions and input. //Author: We need headings to break the article into manageable, logical sections. So, I have inserted headings into what I hope are appropriate places. Feel free to change and move them.//


Understanding turbulence through computation
The computational methods most often used by engineers, even today, are conceptually similar to the methods suggested by the first investigators.In 1897 Joseph Boussinesq suggested replacing the viscosity, a parameter in the Navier-Stokes equation of fluid motion, with a turbulent viscosity that is typically much larger in magnitude and varies with position and time depending on the state of the flow. 3In this and similar methods, the resulting equations are still Navier-Stokes-like, but they differ in the phenomenology of the flow-dependent parameters.Such phenomenology can be refined indefinitely, and ever more experience can be incorporated into it.However, a satisfactory understanding of turbulence still has not been realized, and wind tunnels are still necessary.The reason is probably that these practical, phenomenological methods ignore the many length scales that seem to be an essential, characteristic feature of turbulence.
Turbulent flows are not steady, but fluctuating, and hence some kind of statistical description is necessary.Indeed, the phenomenology mentioned earlier typically involves averages of fluctuating quantities.Perhaps the most striking success of turbulence theory is the prediction by Andrey Kolmogorov of the statistical distribution of energy over the many length scales in turbulent flow, the −5/3 law. 4 By what amounts to simple dimensional analysis, Kolmogorov argued that in strong, fully developed turbulence, the energy distribution function should follow a power law, with more energy in the large-scale motions (large eddies) and less in the small-scale motions, at least over length scales between the large scale, at which the fluid is forced, and the small scale (now called the Kolmogorov length scale), at which viscosity becomes important.Kolmogorov's argument does not require the Navier-Stokes equation at all, so it is not really part of a computational method, but it is a striking insight, nonetheless.The energy in real turbulence does, in fact, seem to be distributed according to the Kolmogorov power law.
Digital computers were applied to the turbulence problem as soon as it was feasible, and new phenomena were discovered soon thereafter.For example, the subject of chaos theory owes much to Edward Lorenz and the famous Lorenz attractor, discovered as a byproduct of trying to model the turbulent atmosphere. 5One reason that chaos became so immediately interesting and popular was the existence of examples such as the Lorenz attractor, involving only a few degrees of freedom and very simple computations.Hence chaos was, in a sense, accessible to everyone, unlike the problem of explicitly simulating turbulence.
Supercomputers of the 1970s finally had sufficient speed and memory for direct numerical simulations of turbulent flow. 6,7These simulations gave researchers a "microscope" to see every detail of turbulence.Initially, the microscope had fairly low resolving power; that is, the grids on which the simulations were performed were fairly crude.With advances in hardware came finer grids and increased resolution.Now these simulations look like real turbulence.For example, the Kolmogorov power law for the energy spectrum can be observed in sufficiently large simulations. 7he supercomputers of the early 1970s are the desktop computers of today.That observation brings us to the substance of this column: most of us have turbulence microscopes sitting on our desks.They are lowresolution microscopes, but they are widely available.Thus the problem of turbulence could be of even more general interest.Increased attention and interest could by itself transform the field.At least this is what we came to believe during a Research Experience for Undergraduates project at Mount Holyoke College in the summer of 1999, where we found that a 200-MHz PC with 32 Mbytes of memory running Matlab could simulate weakly turbulent flow in three dimensions.

Simulating turbulence
The Navier-Stokes equations are Newton's laws of motion applied to a fluid.To derive them one imagines a small volume of the fluid and considers all the forces on it, including most importantly the forces due to the adjoining fluid. 1 Some assumptions must be made, including that there are forces proportional to the rate of strain (viscous forces).The success of the theory in the laminar flow regime for many fluids (Newtonian fluids) is strong evidence that these assumptions are correct.For an incompressible fluid of constant density ρ, the equation for the velocity field u r is ( ) where ν is the kinematic viscosity, P is the pressure, and f r is an external force per unit mass, such as the acceleration due to gravity.
The flow of an incompressible fluid can be described by a divergenceless velocity vector field ( ) , , , u t x y z r , where (x, y, z) is the position.In our simulations positions are defined on a 32 3 spatial grid.(We frequently reduced the grid to 24 3 .)Our aim is to generate turbulent flow, not to model a realistic experimental geometry.We impose periodic boundary conditions on u r , so that there are no solid boundaries.The discretized, periodic flow u r is represented by a finite Fourier series where the sum is over a 32 3 lattice (or 24 3 lattice) in k-space.We also consider the vorticity vector field which is the curl of the velocity vector field.The vorticity has a Fourier representation identical in form to Equation 2, with coefficients Because the flow is divergenceless, we also have 0 , and We introduce a simple force with just two nonzero Fourier components, to do work on the fluid.This force pushes the fluid in the positive x-direction in half the system and in the negative x-direction in the other half, so that the fluid is sheared.In terms of the vorticity ω, the Navier-Stokes equations become The advantage of this form over the more familiar form of Equation 1 is that the pressure P has been eliminated.We express Equation 8 as a system of ordinary differential equations for the time-dependent Fourier coefficients k b r r of the vorticity, where the k c r r in Equation 9are the Fourier coefficients of the nonlinear term, in Equation 8.
The forcing terms k g r r are all zero except for Once we determine the Fourier coefficients k c r r , we have specified the differential equations in Equation 9, and well-known numerical methods can be used to integrate them.(We used a fourth-order Runge-Kutta algorithm.)The method for finding the k c r r proceeds as follows.At each step of the integration of Equation 9, we know k b r r , the Fourier coefficients of w r , and hence by Equation 6, k a r r , the Fourier coefficients of u r .
By inverse fast Fourier transforms, we find u r and w r and hence u w r r . By a fast Fourier transform we find the Fourier coefficients k d r r of this cross product.Then Matlab code for this pseudospectral method can be found at www.mtholyoke.edu/~mpeterso/reu/99/code.
If the viscosity ν is sufficiently large, a random initial flow becomes laminar, because the viscous term  9damps out all the unforced Fourier coefficients.For smaller ν, the nonlinear term k c r r couples energy from the forced coefficients to the other coefficients faster than the viscous term can damp them out, and they do not approach zero.If ν is not too small, a statistically steady state results, in which energy flows from the forced coefficients at small k to the other coefficients, finally to be damped out at the largest k values, where the damping is most effective.These are the turbulent solutions.The flow of energy in k space, accompanied by noisy values for the intermediate Fourier coefficients, is sometimes called the Kolmogorov cascade. 8If ν is chosen to be too small, there is insufficient damping to dissipate the incoming energy, and the numerical solution blows up.
The turbulent solutions in these simple simulations exhibit the characteristic properties of real turbulence.For example, there is anomalous resistance to forcing, as shown in Figure 1.Here the z component of the forced Fourier coefficient b z (0, 1, 0) has been artificially set to zero at time t = 0, and we follow its response in time for a turbulent solution.The dashed line shows the exponential approach to its steady laminar value, the behavior it would have in the absence of turbulence.But because this mode loses energy to other Fourier modes, it never reaches this value, and it levels off noisily at about one-fifth of its laminar steady state value.An experimentalist would say the resistance was anomalously large by a factor of 5. We can see how the energy is distributed over various length scales by summing the energy of the modes in concentric spherical shells in k-space.We slice k-space up into spheres of radius k j about the origin and thickness κ, and find the sums ( ) as shown in Figure 2.This calculation is necessarily rather crude, because we don't have much k-space to work with, just our 32 3 grid.We can't take too many spherical shells, or there wouldn't be enough grid points in each shell to give a good average.However, the distribution of energy looks much the same for any reasonable number of shells and is consistent with a power law, although steeper than the −5/3 Kolmogorov law.The slope in the log-log plot is actually close to −3.This distribution is not unreasonable, because the smallest-length scale must be near the Kolmogorov length in order for dissipation there to establish the statistically steady state.Beyond this length scale (larger k), the energy must fall off faster than any power.We are trapped near this region by hardware limitations, but the distribution still looks like a power law (over an admittedly small range).Our preliminary observations suggest that the qualitative features of turbulence are exhibited in these small simulations and that relatively small simulations might yield physical insight.What such an insight might be is unknown, but it would involve both the geometry and the statistics of the flow.A coordinated search for it is conceivable.The limitations of desktop computers could be partially offset by a distributed approach to such computations.In astronomy, number theory, and cryptology, it has been possible to enlist amateurs on significant projects.Perhaps this approach would also work for turbulence theory.
The simulations that we have described do not correct for aliasing, although aliasing is a problem in the pseudospectral method (the application of finite Fourier transforms to nonlinear problems).Aliasing causes a sufficiently high spatial frequency, generated by the nonlinear term, to be mistaken for a low spatial frequency, because of the discrete sampling imposed by the grid.In our solutions it amounts to a small noise source that is not explicitly included in the model but is an artifact of the method.One can rigorously eliminate aliasing by further truncation of the grid, throwing away the largest k values in each time step, but as Figure 2 shows, these coefficients are small anyway.The review article by Parviz Moin and Krishnan Mahesh is useful on this and other numerical issues. 7See also the book by Claudio Canuto and his colleagues. 9ne might think that the choice of initial conditions would be unimportant and that they could be chosen randomly, but this is not strictly true.If one chooses the k b r r to be random and real, the system is on a rather special, invariant surface in the fluid's phase space, and such flows apparently do not become turbulent!A systematic way to initiate a turbulent flow is to start with what are called ABC flows.These are eigenstates of the curl operator and look typically like A single ABC flow decays exponentially in time, but a superposition of two linearly independent ones becomes turbulent after an initial transient period.

Making sense of the data
Direct numerical simulations of turbulence contain a wealth of information.The challenge is to make sense of it.There are literally thousands of Fourier coefficients, all looking qualitatively like Figure 3.One naturally thinks of using statistical methods.For example, time autocorrelation functions indicate a mixture of randomness and systematic motion.Most correlations decay rapidly, but there is typically an indication of long-term, low-amplitude oscillations as shown in Figure 4.The Kolmogorov picture gives at least a rough picture of what happens in turbulence in Fourier space, but what does it correspond to in real space?Are there just a few important processes, processes that we can visualize?Is there a probabilistic description of this far-from-equilibrium situation that generalizes the Maxwell-Boltzmann distribution?These are completely open questions.Many approaches are possible.
One approach is to try visualization methods.Figure 5 shows a projection of u r onto a typical z = constant slice of space (the z-component is not indicated).The overall mean flow pattern corresponding to the systematic forcing term is not really discernible, because it is complicated by whorls and reconnections that are different in different slices and that change in time.It would be possible to follow the path of a marker advected with the flow and visualize this marked region, analogous to observing dye streaks in a real experiment.One can also construct derived quantities such as the rate-of-strain tensor σ ij .The Fourier coefficients of this quantity are ( ) ( ) ( ) ( ) 1 2 , so that σ ij can be constructed with an inverse Fourier transform.If δ ij is the Kronecker delta, δ ij + σ ij ∆t is the matrix of a quadratic form, defining an ellipsoid at each point of our grid, which is a visualization of the deformation of a small sphere of fluid in the short time ∆t.An eigenvalue of σ ij is positive if the sphere is being stretched along the corresponding principal axis of the ellipsoid, and it is negative if the sphere is being contracted.Because our fluid is incompressible, the trace of σ ij is zero, and both signs typically occur; that is, the fluid stretches along at least one direction and contracts along at least one direction.
The orientation of the vorticity vector w r with respect to the maximal stretching direction has been the subject of much investigation.An early suggestion that vorticity would tend to align with the direction of maximal stretching 10 was shown in direct numerical simulations in the late 1980s to be not generally true, 11 and in fact in certain situations there is a clear tendency for vorticity to align with the eigenvector of the intermediate eigenvalue of the rate of strain. 12This eigenvalue could have either sign, and suggests that what is going on is more subtle than simple, unambiguous stretching of vorticity.
A scalar quantity related to the local rate of stretching of vorticity is We find I to be concentrated, at any given moment, in a few localized regions that show up as "hot spots" (regions where I is large) in slices through the fluid.Another scalar quantity that is nonuniformly distributed in a qualitatively similar way, with local concentrations, is the helicity, u w ¼ r r .Could quantities like these, visible in the turbulence microscope, be keys to a better understanding of turbulence?

Order out of chaos
Over the last few decades, questions that originated in the turbulence problem have given rise to chaos theory, a new paradigm.The Lorenz attractor, distilled out of a system of equations for meteorology (that is, turbulent flow of the atmosphere), was perhaps the first widely discussed example of a simple chaotic system. 5Now it is possible to turn these still relatively new methods back upon the problem that inspired them, turbulence itself.
One simple method for discovering patterns in a chaotic time series {x(t j )} is reconstruction. 13The triplets, (x(t j ), x(t j+k ), x(t j+2k )), are plotted using three axes for some suitable delay k.As an example, a periodic function would generate points that fall on a single closed curve in 3-space.Chaotic functions, although not periodic, might still generate recognizable patterns in space, not so different from closed curves, and readily distinguishable from uniform, random distributions.Figure 6 shows the data of Figure 3 graphed in this way and viewed from two different "directions" in this abstract space.The points of the time series are not randomly distributed but lie predominantly in a plane.Within the plane the curve loops systematically in a way that is not periodic but does not look merely random, either.Such patterns are candidates for one of the basic objects of chaotic dynamics, a "strange attractor."Figure 6.The time series in Figure 3 is reconstructed as described in the article and viewed from two different angles.The points seem to lie approximately in a plane, which is viewed edge-on in the lower graph.The pattern is viewed face-on in the upper graph and shows a repetitive but not periodic behavior, typical of a strange attractor.
A lower bound on the dimension of a strange attractor is the number of positive Lyapunov exponents. 14he existence of at least one positive Lyapunov exponent is the standard criterion for chaos.We estimated the largest Lyapunov exponent 15 for one of our flows and found that the largest Lyapunov exponent is convincingly positive, and hence that our turbulent flows are chaotic.
In this computation we regard a flow at any moment as given by its Fourier coefficients, and hence by a point in the 24 3 -dimensional space of these coefficients.Over time, the flow traces out a trajectory P in this space.We give the space the usual Euclidean norm | |.That is, if {b j } and {B j } are the Fourier coefficients of two flows P b and P B , we say the flows are separated by a distance ( ) By forward integration of Equation 9 from some initial point P(0), we find the trajectory P. At the same time, we start a second trajectory at Q 0 (0), with P(0) − Q 0 (0) = ⑀.(16) After one integration time step, we compute the new distance between the two trajectories: Here λ 1 is the factor by which the initial distance grew, and its logarithm is already a (very crude) estimate for the largest Lyapunov exponent.Because the second trajectory Q 0 might get too far away from P, we start a new trajectory Q 1 by moving the point Q 0 (∆t) toward P(∆t) until it is again just the distance ⑀; that is, we have Then we repeat and find λ 2 from |Q 1 (2∆t and so forth.The largest Lyapunov exponent λ is then ( ) where N is the number of iterations.
Our intuitive picture of a strange attractor suggests that λ should be insensitive to the choices of ⑀, ∆t, P(0), and Q 0 (0), but in reality it is not even clear that the limit exists.Two attempts to calculate the largest Lyapunov exponent are shown in Figure 7.Each calculation seems to be approaching a limit, but they do not appear to be the same limit.In both cases ⑀ = 0.1, but they differ in the choice of P(0), Q 0 (0), and ∆t.The result in both cases is that λ N is convincingly positive.This calculation is an example of one where many desktop computers working separately could perhaps best clarify the pattern, because the computations are time intensive.(Figure 7 represents about two weeks of computation.)Turbulence theory has a well-earned reputation for being difficult, but, paradoxically, it is so difficult that it is actually accessible.Precisely because it is difficult to make progress analytically, an increasingly important part of turbulence theory is computational.The peculiarities of real turbulence suggest that new mathematics and new physics are waiting to be discovered.The power of even modest computers now opens the door to this avenue of investigation.
//Author: Our style is to include the issue number and the month for periodical listings.We'd appreciate this information if you have it on hand.//

Suggestions for Further Study
Chaos theory has introduced the notion of the Lyapunov exponent.In the main article we described an algorithm to compute the largest Lyapunov exponent for our turbulent flows, but what we found was puzzling.It would be useful to know what the "largest Lyapunov exponent" really means for these flows.This

Figure 1 .
Figure 1.The temporal response of the z component of the forced Fourier coefficient

Figure 2 .
Figure 2. A log-log plot of the spectral energy density Φ(k) versus k at a typical moment of the simulation, obtained by summing the energy contained in spherical shells in k-space.The entire accessible range of k was partitioned into 8 (+), 9 (o), and 10 (*) spheres.The three choices gave qualitatively similar results.A line of slope −5/3, the Kolmogorov prediction, is also shown.

Figure 4 .
Figure 4.The autocorrelation function of the Fourier coefficient shown in Figure 3. Initially, the correlation falls off steeply with time, but instead of approaching zero, it oscillates noisily, indicating small correlations over long times.

Figure 5 .
Figure 5. Projection of a typical turbulent velocity field onto a z = constant plane.The z component of the field is not indicated.

Figure 7 .
Figure 7. Two determinations of the largest Lyapunov exponent of one of the turbulent flows.The quantity λ N on the right side of Equation 20 before the limit N → ∞ is taken is shown as a function of t = N∆t.In both cases, ⑀ = 0.1.For ∆t = 0.005, λ N Ϸ 2, and for ∆t = 0.02, λ N Ϸ 13.The two trajectories analyzed here had different starting points.