Intermittency at Fine Scales and Complex Singularities of Turbulent Couette Flow
IIntermittency at Fine Scales and Complex Singularities ofTurbulent Couette Flow
Andre Souza and Divakar Viswanath
Department of Mathematics, University of Michigan (sandre/[email protected]).
Abstract
Fine scales of turbulent velocity fields, beyond the inertial range and well into thedissipative range, are highly intermittent. It has been hypothesized that complex planesingularities are the principal mechanism behind fine scale intermittency. In this article,we view the velocity field of a turbulent flow as an analytic function of time. Althoughthe function is only available for real values of time, we present a numerical technique toanalytically continue the function to complex values of time, and with sufficient fidelity tolocate and visualize the singularity closest to the real axis. Using this technique, we demon-strate a robust connection between temporal intermittency and the location of singularitiesin the complex plane.
Intermittency in turbulent flows refers to sharp and spatially isolated peaks of energy containedin large wave numbers, beyond the inertial range and well into the dissipative range. In thewords of Batchelor and Townsend [2], who discovered this phenomenon, “the energy associatedwith large wave-numbers is very unevenly distributed in space” and “there appear to be isolatedregions in which the large wave-numbers are ’activated’, separated by regions of comparativequiescence.” This fine scale intermittency is a fundamental feature of turbulent flows [1, 5, 17].Frisch and Morf [7] hypothesized that intermittency is a manifestation of singularities in thecomplex plane. Thus the activation mechanism is believed to be the occurrence of singularitiesin the complex plane, with the intermittent peak occurring right below the singularity closest tothe real axis. Complex singularities are intrinsically easier to study with a single independentvariable. Therefore Frisch and Morf shifted the emphasis to temporal intermittency, with timeas the single independent variable with the entire velocity field viewed as a function of time.The connection between intermittency and complex singularities has been demonstrated ina variety of systems such as Langevin’s and Burger’s equations [7, 15]. However, to the bestof the authors’ knowledge, such a connection is yet to be demonstrated for the incompressibleNavier-Stokes equations in the turbulent regime. In this article, we develop a numerical tech-nique that is applicable to trajectories of large scale PDE. Using that technique, we demonstratea robust connection between temporal intermittency and complex singularities.Here at the outset, we reprise a beautiful argument of Kraichnan [12]. The argument isheuristic but gives powerful intuition regarding the phenomenon of fine scale intermittency,1 a r X i v : . [ phy s i c s . f l u - dyn ] O c t artially explaining many of the plots found later in this article and showing why the manifes-tation of intermittency is clearest in the dissipative range. The argument is as follows. Supposethat the velocity field of the turbulent flow is infinitely differentiable with each derivative beingsquare integrable. Direct numerical simulation of turbulence from the incompressible Navier-Stokes equations implicitly assumes the velocity field to be analytic, and the great successof direct numerical simulation is certainly evidence for the validity of that assumption. Theassumption of smoothness implies that energy should fall off faster than algebraically (andexponentially fast if the velocity field is analytic) with increase in wave-number. The super-algebraic and possibly exponential decay of energy with wave-number means that a change inwave-number “by a few percent” implies an “enormous” change in the energy associated withthe wave-number. If the energy spectra is governed by a length scale that varies gently fromregion to region, this gentle variation is enormously amplified and shows up as intermittencyin the high wave-numbers. In essence, the hypothesis of Frisch and Morf [7], investigated andvalidated in this article, is that the gentle variation in the large which gets amplified intointermittency in fine scales is due to singularities in the complex plane.Energy dissipation per unit mass, denoted (cid:15) , is an important parameter in turbulence.Indeed, it is the only parameter in the Kolmogorov theory of the inertial range [1]. Fine scaleintermittency may be detected experimentally through an analysis of the variation in (cid:15) [13, 14].Thus physical space analysis of turbulent velocity fields obtained from numerical simulationsuffices for demonstration of spatial intermittency as well as for isolating structures associatedwith that phenomenon [16].The tack we take in this article is different, and possibly complementary. We study finescale intermittency through analytic continuation into the complex t -plane. An advantage isa direct connection to the hypothesized mechanism. Continuation into the complex t -planecould be complementary to physical space analysis as we discuss in the concluding section.Analytic continuation into the complex plane is infamous for its numerical difficulty. Thetechnique of Padé approximation, and strategies for better numerical stability, have been setforth by Weideman [25], Gonnet, Pachón, and Trefethen [8] and Webb [22]. A crucial ingredientis the use of complex phase plots championed by Wegert and others [23, 24], and whosepertinence to the computation of complex singularities was communicated to us in person byHrothgar and Trefethen. These phase plots play a central role in our technique and we rely onthem to locate singularities in the complex plane. The phase plots proved to be a more viableapproach to stabilize numerical computations than strategies for deflating tiny and repeatedsingular values developed in [8].All numerical simulations are carried out with constant time steps for reasons explainedin section 2. The numerical method for locating complex singularities described in section 2interpolates to Chebyshev points of the second kind from equi-spaced data. The Chebyshevseries are filtered before computing Padé approximants. Many of the ingredients in the nu-merical method we describe have occurred in earlier work in some form. Our contribution isto synthesize a numerical method that works robustly for turbulent flows.Somewhat counter-intuitively, the high dimensionality of discretizations of turbulent veloc-ity fields appears to be a help rather than a hindrance. The numerical method of section 2takes advantage of intermittency in high wave-numbers in locating the complex singularities.In section 3, we compute the complex singularities of six periodic or relative periodicsolutions of plane Couette flow. The plane Couette flow set-up we use is the minimal flow unit210]. The minimal flow unit was derived by constraining turbulent Couette flow into a small boxat the low Reynolds number of Re = 400. Since the Reynolds number is low, the dissipativerange is reached easily. Another advantage of the minimal flow unit is its connection to theself-sustaining process developed by Waleffe [21]. While inertial range turbulence is statistical,turbulence is also characterized by coherent structures. The self-sustaining process emphasizesthe dynamical aspects of fine-scale turbulence with connections to coherent structures.The six periodic and relative periodic orbits, whose complex singularities are computed, arefrom earlier work [19] . For each of these orbits, we demonstrate a clear connection betweenintermittency and the location of complex singularities. There are no profound advantagesto using periodic orbits. One could use any turbulent trajectory. But as a practical matterthere are two reasons for opting for these solutions of plane Couette flow. Firstly, the orbits arereported in [19] with estimates for their precision so that they may be reproduced. Secondly, theearlier study of these orbits has identified different regimes enabling a preliminary investigationof the connection between observed turbulence and the location of complex singularities.The concluding section 4 discusses further connections of the work presented in this article.One of these is to the spatial analysis of fine scale intermittency, already alluded to above.Another is to blow-ups in physical space. The numerical method used for computing singularities will be described using the Lorenzexample. The Lorenz equations ˙ x = 10( y − x ) , ˙ y = 28 x − y − xz, ˙ z = − z/ xy admitsingular solutions that may be expanded in convergent psi-series as follows [20]: x ( t ) = P − ( η ) t − t + P ( η ) + P ( η )( t − t ) + P ( η )( t − t ) + · · · y ( t ) = Q − ( η )( t − t ) + Q − ( η ) t − t + Q ( η ) + Q ( η )( t − t ) + Q ( η )( t − t ) + · · · z ( t ) = R − ( η )( t − t ) + R − ( η ) t − t + R ( η ) + R ( η )( t − t ) + R ( η )( t − t ) + · · · (2.1)Here η = log( t − t ) with the branch cut chosen appropriately, and P i , Q i , R i are polynomials in η . These polynomials are given explicitly in earlier work [20]. Here we note that Q − = − i/ R − = 1 / P − = Q − = 2 i , and R − = 17 /
9, so that the dominant part of the singularity isa pole for x ( t ) and a double pole for y ( t ) and z ( t ). The first terms in the psi-series which arenot constants are Q − = − i − i ( η + C ) , R − = 138554 − η + C ) , where C is an undetermined constant. There is another undetermined constant D in thepsi-series which occurs for the first time in Q .Not all the singularities of Lorenz are proved to be psi-series of this form. However, a fewsingularities were examined numerically in [20] and were all found to be psi-series of the formgiven above. In addition, the locations of the singularities were determined with 10 digitsof precision, which makes the Lorenz example useful for validating the numerical methoddescribed here. 3he Lorenz example, being of only 3 dimensions, can be tackled using a number of tech-niques. It is amenable to extended precision computations. Orbits of the Lorenz system werecomputed with 500 digits of accuracy in [20]. For the incompressible Navier-Stokes equations inthe fully turbulent regime, simply integrating the trajectories is among the most complex anddifficult computations. The room for deploying techniques that find complex plane singularitiesis much more limited.In this section, we refine Padé based techniques [8, 25] and derive a numerical methodthat is reliable when applied to turbulent signals. The numerical method is described in threestages. The first stage is an analysis in the purely real time domain. The second stage is thecomputation of a Padé approximant. The third stage is a discussion of the Wegert–style phaseplots of the Padé approximant. Complex singularities are located using phase plots, and wecompare the locations obtained through the phase plots to the far more accurate computationsof [20]. To describe the numerical method, we look at the z -coordinate of the AABB orbit of theLorenz system. For the nomenclature of Lorenz orbits, see [18]. The
AABB orbit has periodequal to 3 . ... The Lorenz orbits can be computed using a special method described in [18]. Turbulencetrajectories are obtained using ordinary direct numerical simulation, and to obtain a reliablecomparison the Lorenz orbits too were computed using numerical integration. There is oneimportant point pertinent to direct numerical simulation, however. The time-steps must beconstant.If a Hamiltonian system is integrated using constant time-steps and a symplectic discretiza-tion, it has been proved that the numerical trajectory (ignoring rounding errors) approximatesthe trajectory of a perturbed Hamiltonian with error that is exponentially small in the time-step for an interval of time that is exponentially long in the time-step [9]. The arguments of[9] can be made to apply to non-Hamiltonian ordinary differential equations without majormodifications. Thus if the time-step is constant, there is reason to think that the effect ofthe discretization error is to introduce a smooth and analytic perturbation of the underlyingdifferential equation. On the other hand, adaptive time-stepping strategies are non-smooth,and even though they allow for longer time-steps, they destroy the analytic structure of theunderlying differential equation.Once a signal u ( i ∆ t ) is obtained for integer values of i , the next step is to pick a timeinterval and compute a Chebyshev series. Suppose that the time interval is [ a, b ] with a = i a ∆ t and b = i b ∆ t . To convert, the signal u ( t ), a ≤ t ≤ b , into a Chebyshev series, we firstinterpolate the signal at Chebyshev points (of the second kind) cos( jπ/n ), j = 0 , , . . . , n ,shifted for the interval [ − ,
1] to the interval [ a, b ]. If t j is a shifted Chebyshev point in [ a, b ]and t j ∈ [ i ∆ t, ( i + 1)∆ t ], the value of u ( t j ) is computed through polynomial interpolation usingthe nodes t = ( i − k ) ∆ t, ( i − k + 1)∆ t, . . . , ( i + k + 1)∆ t .The polynomial interpolation is from equi-spaced data. However, there is no need to fearthe Runge phenomenon as t j is near the center of the domain of interpolation. In addition, k is small, with k = 4 in almost all reported computations ( k = 4 implies at least 8-th orderaccuracy). The interpolation was carried using the barycentric Lagrange formula, which isknown to have excellent numerical stability [3]. Interpolated Chebyshev data was moved back4 .0 0.2 0.4 0.6 0.8 1.0 1.210152025303540450 10 20 30 40 50 60 7010 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 -9 -8 -7 -6 -5 -4 -3 -2 -1 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 Figure 2.1: Plots of the z -coordinate of the AABB
Lorenz orbit. The top plots show z as afunction of time, with the dots being interpolants at Chebyshev points. The bottom plots are | c n | , absolute values of coefficients of the Chebyshev series, of the top plots graphed against n .to the original equi-spaced grid to assess the quality of interpolation, and to ensure that noinformation was lost during interpolation.Once the signal u ( t ) is available at Chebyshev points shifted to the interval [ a, b ], it isconverted to a Chebyshev series using the discrete cosine transform. If T k ( x ) = (cid:16) z k + z k (cid:17) denotes the Chebyshev polynomial of degree n , with z determined from x = (cid:16) z + z (cid:17) , and ˜ T k is T k shifted from the interval [ − ,
1] to the interval [ a, b ], the Chebyshev series obtained fromthe discrete cosine transform is of the form u ( t ) = c + c ˜ T ( t ) + c ˜ T + · · · + c n ˜ T n and is valid for a ≤ t ≤ b . This Chebyshev series will be converted into a Padé approximantlater. It is important to do a preliminary analysis of the coefficients c i before computing thePadé approximant. Much of the analysis is visual as we will now explain.Figure 2.1 shows plots for the Lorenz orbit AABB . The period of the orbit has beensplit into three equal segments. It is easily noticeable that the Chebyshev series of the middlesegment is cleaner with a more pronounced pattern. The Chebyshev coefficients decreaseexponentially all three segments, but in the middle segment we see spikes that protrude downwith a clear rhythm. Isolating a Chebyshev series with this kind of pattern is crucial. ThePadé approximant relies on this kind of a pattern to locate the complex singularity. If thepattern is not clear to begin with the Padé approximant will not work too well.For some of the Lorenz orbits and even turbulent orbits of plane Couette flow, we are ableto find Chebyshev series that exhibit patterns that are even cleaner and more pronounced thanthe middle segment in Figure 2.1. For others, the pattern is not as clean. In addition to thepattern being clean and well-established, it is advantageous if the pattern starts early.Once a Chebyshev series with a clean pattern has been isolated, the pattern can be improvedfurther using filtering as shown in Figure 2.2. Here the Chebyshev series of the middle segmentwas filtered by setting to zero the first 10 coefficients, since they are not part of a pattern, aswell as the last few which were deemed to be not sufficiently accurate. The filtered signal onthe left side of Figure 2.2 almost gives away the real part of the location of the singularity as < ( t ) ≈ .
6. Filtering makes the intermittency more pronounced.5 .0 1.2 1.4 1.6 1.8 2.0 2.2t2.01.51.00.50.00.51.01.52.02.5 z
10 15 20 25 30 35 40 45 50 55n10 -7 -6 -5 -4 -3 -2 -1 c o e ff Figure 2.2: Filtered signal corresponding to the middle segment of Figure 2.1 and a linear fitof its Chebyshev coefficients.Identification of the role of filtering in accentuating intermittency is due to Frisch and co-workers [7, 15]. As noted by those authors, the roots of the idea are found in the theory ofcomplex variables as developed in the late 1800s and early 1900s. In that body of work, theasymptotics of the coefficients of a power series are connected to the form and location of thesingularities. The further out one goes into the series, the clearer the asymptotics will be. Innumerical work, however, one cannot go too far because of accuracy issues.Alternatively, the role of filtering may be understood as follows. Suppose that f ( t ) isanalytic except for singularities at t = a ± ib , with b = 0, in the complex plane. Let t ∈ R .Consider the expansion of f ( t ) in powers of ( t − t ) with t varying along the real line. Theradius of convergence is least at t = a and the coefficient of ( t − t ) n spikes at t = a , for large n . The increase in steepness of the spike at t = a (in a relative sense) as n increases is thesubstance of the argument for the connection between intermittency and complex singularities.In this setting, it is quite clear that we can add a polynomial or any entire function to f ( t )without moving the singularity. Such an addition can destroy intermittency in the coefficientsof ( t − t ) n for small n , but the intermittency will show up for larger n . Thus intermittencymay be identified more clearly in the signal f ( t ) if the coefficients of ( t − t ) n for certain small n are filtered out. In the numerical method we are describing, the range of filtering is determinedby examining the coefficients of the Chebyshev series.The right hand plot in Figure 2.2 shows a linear fit of the filtered Chebyshev coefficients.The linear fit is obtained as follows. An index i is selected if and only if it is in the range thatsurvives filtering (10 < i <
55 in Figure 2.2) and | c i | > | c j | for j > i . The linear fit of log | c i | vs i is computed using only the selected indices i . The linear fit approximates the outer envelopeof log | c i | vs i as shown in Figure 2.2. If the slope of the fit is ρ , < ( t ) = 12 (cid:18) ρ + 1 ρ (cid:19) cos θ, = ( t ) = 12 (cid:18) ρ − ρ (cid:19) sin θ (2.2)is an ellipse in the complex plane parametrized by θ . Standard results in complex functiontheory imply that there must be a singularity on this ellipse. The Padé approximant is farmore reliable inside the ellipse than outside it, as we will momentarily see.6
10 20 30 40 50 60n10 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 c o e ff Figure 2.3: Comparison of filtered coefficients of Figure 2.2 with coefficients of the Padé ap-proximation. The Padé approximations are dotted.
Given a Chebyshev series, its Padé approximant n X i =0 c i T i ( t ) ≈ a + a T ( t ) + · · · + a p T p ( t ) b + b T + · · · + b q T q ( t )is computed using( c + c T + · · · + c n T n ) ( b + b T + · · · + b q T q ) − ( a + a T + · · · + a p T p ) ≈ . If Chebyshev polynomials are multiplied using the identity T i T j = (cid:16) T i + j + T | i − j | (cid:17) and co-efficients are equated to zero, we get a set of n + 1 linear relations between the unknowncoefficients a i and b i .We follow Gonnet, Pachón, and Trefethen [8] in using the singular value decomposition ofthe last ( n − p ) of the equations to solve for the b i with the normalization P | b i | = 1. Thesame authors recommend strategies to reduce the degree of the dominator if the singular valuesare too small or if the smallest singular values are repeated. We have experimented with thesestrategies, but rely on different criteria to obtain robust Padé approximants. Removing smallsingular values was problematic and threw away too much information. In fact, in most of ourcomputations we prefer to have denominators of large degree with the degree being nearly 200in some cases.Figure 2.3 compares the filtered coefficients and the Chebyshev coefficients of the Padéapproximation. The Padé approximation was obtained using p = 15 and q = 25. Since thefirst 10 coefficients are filtered out, the first 10 coefficients in the numerator are 0. Mostof the fitting is taking place in the denominator. Close observation shows that the Padéapproximation gets the pattern near the tail slightly differently. However, the quality of theapproximation is excellent through the range of the Chebyshev coefficients. Obtaining a highquality approximation of this sort is essential. 7 .0 1.2 1.4 1.6 1.8 2.0 2.20.40.30.20.10.00.10.20.30.4 1.0 1.2 1.4 1.6 1.8 2.0 2.20.40.30.20.10.00.10.20.30.4 Figure 2.4: Phase plots of Padé approximants of the z and x coordinates, respectively, of theLorenz orbit AABB corresponding to the middle segment in Figure 2.1. The colors code forthe phase or argument at a point in the complex t -plane.If two analytic functions agree exactly for t ∈ [ a, b ], they will agree everywhere in thecomplex plane. If one of those two functions is a polynomial, it is incapable of picking upsingularities of the other function, even if the agreement over t ∈ [ a, b ] is very close, becausepolynomials are entire functions. The hope with Padé approximations is that they will pickup the singularity in the complex plane because they are rational in t . For that to be true, weneed close agreement over t ∈ [ a, b ] at a minimum. In fact, such close agreement also appearsto be sufficient if the original signal is properly filtered. As mentioned above, the left side plot in Figure 2.2 gives a good indication of < ( t ) if t is thelocation of the singularity. The imaginary part = ( t ) may be obtained by plotting the Padéapproximant for complex values of t . Phase plots [23, 24] of a complex function f ( z ) depictthe function by indicating the phase arg f ( z ) using a color code. Such plots reveal a wealth ofinformation about poles and zeros. According to Wegert and Semmler [24], the precise originof the idea of phase plots is difficult to determine.Figure 2.4 shows phase plots for the z (left side) and x (right side) coordinates of the middlesegment of the Lorenz orbit AABB . Each of the 11 zeros of the filtered signal shown in Figure2.2 may be easily located in the phase plot on the left. At a zero the colors red-blue-greenrotate in a counter-clockwise sense (opposite of Wegert’s convention). The ellipses shown asthick lines correspond to (2.2).The phase plot also picks up a singularity near the boundary of the ellipse and with < ( t ) ≈ . z ( t ), the Padé approximantappears to be attempting to pick up a double pole. The colors outside the ellipse are notas reliable as those inside. Yet one may see that the Padé approximant is attempting twoclockwise rotations of the colors red-blue-green. In the right hand plot, corresponding to x ( t ),there is clearly a single pole at the same location. The double pole in the z ( t ) phase plot is8rbit x y z accurate AB .
175 0 .
176 0 .
172 0 . AAB .
170 0 .
167 0 .
169 0 . AAAB .
162 0 .
158 0 .
150 0 . AABB .
170 0 .
170 0 .
156 0 . t -axis for four Lorenz orbits.The x, y, z columns are from singularity locations determined using the x, y, z coordinates,respectively. The final column is from extended precision computations [20].Orbit Re τ L M N dx + dy + dz + P1 28 32 64 48 4 . . . . . . . . . . . . . . . . . . z ( t ).Similarly, the phase plot for x ( t ) is picking a single pole because the dominant part in thepsi-series is a single pole. In both cases, the artifacts of Padé approximation are well outsidethe ellipse.Table 2.1 furthers validates the numerical method by comparing the complex singularitylocation determined using x , y , and z coordinates with a more accurate computation fromearlier work. It appears that the singularity location is determined correctly with about twodigits of precision. Since the method for determining the singularity is graphical, and basedon mouse-clicks, two digits of precision is about as much as we expect.We have avoided finding roots of the polynomial in the denominator of the Padé approxi-mant. Polynomial root finding is often numerically unsound. For the Lorenz orbits, denomina-tor degrees range from 25 to 40. For the Couette orbits that we next turn to, the denominatordegrees can be larger than 100 and polynomial root-finding is infeasible. The phase plots aremore reliable and informative. The velocity field of plane Couette flow is represented as u ( x, y, z, t ) = L/ X ‘ = − L/ N/ X n = − N/ M X m =0 u l,m,n ( t ) exp (cid:18) i‘x Λ x + inz Λ z (cid:19) T m ( y ) . Thus the velocity field is periodic in x and z with periods equal to 2 π Λ x and 2 π Λ z . Forthe minimal flow unit, these periods are 5 . .. and 3 . ... , respectively [10]. As before,9 m are Chebyshev polynomials. At the walls y = ±
1, the velocity field satisfies the boundarycondition u = ( ± , , ∂ u ∂t + ( u . ∇ ) u = −∇ p + 1 Re u is Re = 400. The pressure p is determined by the incompressibility constraint ∇ . u = 0.Table 3.1 shows resolution data for six periodic or relative periodic solutions of the minimalflow unit. A more detailed description of these orbits may be found in [19]. The orbits arelabeled P1 through P6. P1 was first computed by Kawahara and Kida [11].The frictional Reynolds number Re τ of all six orbits is low, although they exhibit somefeatures of turbulence such as a sharp peak in turbulence energy production in the near wallregion. The orbits capture other aspects of buffer layer dynamics as well [19]. The purpose ofthe minimal flow unit is to capture certain essential features of turbulence in as small a systemas possible. The fact that the Reynolds number is low, makes it easier to reach into the highwave-number region where intermittency manifests itself most clearly.Table 3.1 gives the spacing between grid points in the x, y, z direction in frictional units.For the non-uniform Chebyshev grid in the y -direction, the reported number is the maximumspacing. There is no question that all six orbits are exceedingly well-resolved. Good resolutionis essential for reaching wave-numbers where intermittency is pronounced.The Lorenz system is only 3 dimensional. For the minimal flow unit, the dimensionality ofthe spatial discretizations we use is more than 10 . Since the numerical method of the previoussection is visual, only a small fraction of these modes can be examined. We examine about 50modes for each orbit. These modes are indexed by k . Let δ‘ = L/ δm = M/
32, and δn = N/
64. For 0 ≤ k mod 16 < k = 15, these modes are ‘, m, n ≈ jδ‘, jδm, jδn ,where j = k mod 16. For 6 ≤ k mod 16 <
12, we use ‘ = m = n = j −
5, where j = k mod 16 as before. For 12 ≤ k mod 16 <
15, we use ‘, m, n = L/ j − , M/ j − , N/ j − j = k mod 16 as before.If u = ( u, v, w ), we extract the stream-wise u mode if 0 ≤ k <
16, the wall-normal v modeif 16 ≤ k <
32, and the span-wise w mode if 32 ≤ k <
48. Most of these modes are complexvalued, and we retain only the real part in the signal that is extracted.The selection of modes for 0 ≤ k <
48 is admittedly a little arbitrary. The aim is to look ata few low wave-numbers, a few in the middle, and a few higher wave-numbers. The choice ofwave-numbers proved adequate for demonstrating the connection between intermittency andcomplex singularities in all six orbits.Table 3.2 shows the parameters used for computing the complex singularities. The mode k and the filtering interval are chosen so that the coefficients of Chebyshev series in time exhibita clean pattern. Plots of the filtered series are shown in Figure 3.1 for P1, P2, P3, and P5. Ineach case, a clear pattern is discernible, and the Padé approximants reproduce the coefficientswith excellent accuracy.The plots in Figure 3.1 pick modes corresponding to values of k listed in Table 3.2. Wetried to pick the k for which the pattern in the Chebyshev series was the cleanest. Althoughthere are other choices of k that are nearly as good, a pattern shows up cleanly only for a fewvalues of k .The need to select k carefully may be explained as follows. In complex function theory, f ( z ) is usually a complex valued function of a complex variable. However, much of complexfunction theory goes through with little change even if f ( z ) takes values in a Hilbert space or10rbit k n fst n lst p q d P1 19 12 65 18 45 12 . . . . . . k column gives the mode for which theChebyshev series in time had the best pattern. The interpretation of k is found in the text.All Chebyshev modes c i outside of ( n fst , n lst ] were filtered out. The p and q columns givethe degrees of the Padé numerator and denominator, respectively. Notice that most of thecoefficients in the numerator are nearly zero because of filtering. The last column gives theestimated distance of the singularity from the real t -axis. -23 -21 -19 -17 -15 -13 -11 -9 -7 c o e ff -20 -19 -18 -17 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 c o e ff -20 -19 -18 -17 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 c o e ff -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 c o e ff Figure 3.1: The filtered Chebyshev series for P1, P2, P3, and P5 (clockwise from top-left).The y -axis is the absolute value of the Chebyshev coefficient of T n ( t ). The circular markersare approximations computed using the Padé approximant.11 Banach space, as long as z is a complex variable. For instance, Cauchy’s theorem continuesto be true.We may think of u ( t ), the velocity field, as an analytic function of t which takes valuesin some Hilbert space, although the velocity field is computed only for real values of t . Thesingularity at a point t in the complex t -plane may be thought of as a series of some type in( t − t ), something like the psi-series in (2.1), with coefficients in a Hilbert space.When we pick a mode corresponding to some k and record the variation of that modein time, we are essentially taking an inner product of the mode and the expansion of thesingularity. The Padé approximant tries to pick up the dominant part of the singularity. Tworandom vectors in high dimensional space are nearly orthogonal. Ordinarily the dominant partof the singularity does not project well when an inner-product is taken. Thus the need to pick k carefully. t -plane The connection between complex singularities and intermittency for solutions of the minimalflow unit is demonstrated in Figure 3.2. For each of the four orbits shown, the top plot isan intermittent high wave-number signal that was not filtered. The mode that was chosencorresponds to k with 12 ≤ k mod 16 <
16, with the interpretation of k given above. Almostany k in this range can be picked.All the phase plots, except for P1, are clean inside the ellipse. For P4 and P5, thereis a singularity near the boundary of the ellipse, corresponding to the intermittency, but nospurious singularities are found inside the ellipse. The intermittency in P2 is not as pronouncedand there may be more one singularity near the border of the ellipse. The filtered signalscorresponding to the phase plots are not shown. Yet the oscillation of these signals with t forreal t is clear from the phase plots. When interpreting the phase plots, it must be borne inmind that a counter-clockwise cycle of red-blue-green is a zero, while a clockwise cycle is apole .Following [7, 12], the reason that intermittency manifests itself in the high wave-numbers,without filtering, but is more pronounced only after filtering in the lower wave-numbers maybe explained as follows. The energy in the spatial modes of the velocity field decreases expo-nentially with wave-number in the dissipative range. At a given point in time, with t real, thisrate is related in some way to the distance to the nearest singularity in the complex t -plane.The rate appears to be highest right below the singularity. A small change in the rate causes amuch bigger change in the relative magnitude of the energy in the high wave-numbers, makingintermittency quite pronounced as evident from Figure 3.2. We do not claim much accuracy for the distances of the singularities from the real t -axisreported in Table 3.2. The errors may be as high as 20%. Nevertheless, it is very clear thatP1’s singularity is much farther way than that of P2 through P6.Within the confines of the minimal flow unit, P2 through P6 are all in the turbulent regime.However, P1 is intermediate between laminar flow and turbulent flow [11, 19]. The laminar This is the opposite of Wegert’s convention. Figure 3.2: An intermittent high-wavenumber signal as a function of time and a phase plot ofthe Padé approximant after filtering for some other mode (the mode corresponding to k listedin Table 3.2) for the orbits P1, P2, P4, and P6 (clockwise from top left).13olution of course is constant in time and therefore an entire function. P1’s singularity seemsto be much farther out because it is intermediate between laminar flow and turbulent flow.Another question we may ask is about the type of the singularities. Are they poles ofsome order, or are they approximable by a pole of some order? For Lorenz, we know thatthe psi-series are convergent singular solutions. Although not all singularities are proved to beof that form, a few singularities have been analyzed in extended precision and verified to bepsi-series [20]. The dominant part in the psi-series is either a pole or a double pole. Thus thesuccess of Padé is not a big surprise.In the case of the incompressible Navier-Stokes equations, the analytic form of singularsolutions is not known. The Padé approximants use poles to approximate all kinds of singu-larities and do not give direct information about the form of the singularities. Yet the claritywith which Padé approximants locate the singularities, as shown by the phase plots of Figure3.2, makes us think that the form of the singularities is unlikely to be much more complexthan for Lorenz.Another related question is about the number of singularities. Each of the orbits P1 throughP6 seems to be governed by a single dominant singularity, or a tight cluster of singularitiesclosest to the real axis in the case of P2. For a long trajectory and in the complex t -planeas a whole, it appears reasonable to suspect that singularities are numerous and occur withincreasing density as the Reynolds number increases. Computations described in the previous section, even though limited to the minimal flow unitof plane Couette turbulence, offer definite evidence that intermittency in fine scales is due tosingularities in the complex plane as suspected by Frisch and others [5, 7].All our computations of complex singularities are in the time domain. However, intermit-tency also has a significant spatial aspect [17]. The connection between singularities in thetime domain and spatial aspects of intermittency is not entirely clear. It may be possibleto combine spatial modes in a manner that the temporal intermittency becomes particularlypronounced.On the other hand, the velocity field at a fixed value of time may have singularities if thespatial coordinates are allowed to be complex. It is unclear what such spatial singularities looklike and how they relate to temporal singularities.All six orbits discussed in the previous section go through one cycle of the self-sustainingprocess in a single period [19]. The self-sustaining process as developed by Waleffe [21] ischaracterized in the spatial domain in terms of rolls and streaks. Each of the six orbits of theminimal flow unit is governed by a single dominant singularity (and its complex conjugate),or perhaps by a tight cluster of singularities. The possibility of a connection between theself-sustaining process and the location of singularities in the complex t -plane may be worthinvestigating.The finite time blow-up problem for the incompressible Navier-Stokes equations has beenviewed in light of complex singularities [4, 6]. The robustness with which complex singularitiesare located and visualized in Figure 3.2 suggests that computing the distance of the singularitiesfrom the real t -axis as a function of the Reynolds number may shed some light on this topic.However, there is a subtle difference between mathematical and computational investigations14f finite time blow-up.The evidence that direct numerical simulation algorithms solve the Navier-Stokes equa-tions is overwhelming. Yet there is a subtle difference between numerical and mathematicalinvestigations. In mathematical formulations, the velocity field of the Navier-Stokes equationsis allowed to belong to certain function spaces, and these function spaces are much wider thanthe class of velocity fields that are subject to numerical simulation. Any direct numerical sim-ulation code will become unstable if the initial velocity fields are not sufficiently regular. Athigh Reynolds numbers such as Re τ = 1000, establishing a turbulent initial state can pose fargreater difficulties than sustaining turbulence. In fact, at such Reynolds numbers, attemptingto transition to turbulence by adding a small disturbance to laminar flow is nearly impossiblebecause of the large number of transitional instabilities. A turbulent initial condition is estab-lished by slowly raising the Reynolds number. Thus turbulent velocity fields are pre-selectedin a way to induce statistical steadiness in the flow with respect to time integration. We are very thankful to Nick Trefethen, and his students Hrothgar and Marcus Webb, forsharing their work on Padé approximation with us. This research was partially supported byNSF grant DMS-1115277.
References [1] G.K. Batchelor.
The Theory of Homogeneous Turbulence . Cambridge University Press,1953.[2] G.K. Batchelor and A.A. Townsend. The nature of turbulent motion at large wave-numbers.
Proceedings of the Royal Society of London, Series A , 199:238–255, 1949.[3] J.-P. Berrut and L.N. Trefethen. Barycentric Lagrange interpolation.
SIAM Review ,46:501–517, 2004.[4] A. Biswas and C. Foias. On the maximal space analyticity radius for the 3d Navier-Stokesequations and energy cascades.
Annali di Matematica , 193:739–777, 2014.[5] U. Frisch.
Turbulence . Cambridge University Press, 1995.[6] U. Frisch, T. Matsumoto, and J. Bec. Singularities of the Euler flow? not out of the blue!
Journal of Statistical Physics , 113:761–781, 2003.[7] U. Frisch and R. Morf. Intermittency in nonlinear dynamics and singularities at complextimes.
Physical Review A , 23:2673–2705, 1981.[8] P. Gonnet, R. Pachon, and L.N. Trefethen. Robust rational interpolation and least squares.
Electronic Transactions on Numerical Analysis , 38:146–167, 2011.[9] E. Hairer and C. Lubich. The lifespan of backward error analysis for numerical integrators.
Numerische Mathematik , 76:441–462, 1997.1510] J.M. Hamilton, J. Kim, and F. Waleffe. Regeneration mechanisms of near-wall turbulencestructures.
Journal of Fluid Mechanics , 287:317–348, 1995.[11] G. Kawahara and S. Kida. Periodic motion embedded in plane Couette turbulence: re-generation cycle and burst.
Journal of Fluid Mechanics , 449:291–300, 2001.[12] R.H. Kraichnan. Intermittency in the very small scales of turbulence.
Physics of Fluids ,10:2080–2083, 1967.[13] C. Meneveau and K.R. Sreenivasan. Simple multifractal cascade model for fully developedturbulence.
Physical Review Letters , 59:1424–1427, 1987.[14] C. Meneveau and K.R. Sreenivasan. The multifractal nature of turbulent energy dissipa-tion.
Journal of Fluid Mechanics , 224:429–484, 1991.[15] W. Pauls and U. Frisch. A Borel transform method for locating singularities of Taylorand Fourier series.
Journal of Statistical Physics , 127:1095–1119, 2007.[16] J. Schumacher, K.R. Sreenivasan, and P.K. Yeung. Very fine structures in scalar mixing.
Journal of Fluid Mechanics , 531:113–122, 2005.[17] K.R. Sreenivasan. Fractals and multifractals in fluid turbulence.
Annu. Rev. Fluid Mech. ,23:539–600, 1991.[18] D. Viswanath. Symbolic dynamics and periodic orbits of the Lorenz attractor.
Nonlin-earity , 16:1035–1056, 2003.[19] D. Viswanath. Recurrent motions within plane Couette turbulence.
Journal of FluidMechanics , 580:339–358, 2007.[20] D. Viswanath and S. Şahutoğlu. Complex singularities and the Lorenz attractor.
SIAMReview , 52:294–314, 2010.[21] F. Waleffe. On a self-sustaining process in shear flows.
Physics of Fluids , 9:883–900, 1997.[22] M. Webb. Computing complex singularities of differential equations with Chebfun.
SIAMUndergraduate Research Online
Visual Complex Functions: An Introduction with Phase Portraits . Birkhauser,2012.[24] E. Wegert and G. Semmler. Phase plots of complex functions: a journey in illustrations.
Notices of the AMS , 58:768–780, 2010.[25] J.A.C. Weideman. Computing the dynamics of complex singularities of nonlinear PDEs.