On modelling transitional turbulent flows using under-resolved direct numerical simulations: The case of plane Couette flow
aa r X i v : . [ n li n . PS ] O c t On modelling transitional turbulent flows usingunder-resolved direct numerical simulations:The case of plane Couette flow
Paul Manneville and Joran RollandLaboratoire d’Hydrodynamique´Ecole Polytechnique, F-91128 Palaiseau, FranceAccepted on August 24, 2010, to appear in TCFD
Abstract
Direct numerical simulations have proven of inestimable help to our understand-ing of the transition to turbulence in wall-bounded flows. While the dynamics of thetransition from laminar flow to turbulence via localised spots can be investigatedwith reasonable computing resources in domains of limited extent, the study of thedecay of turbulence in conditions approaching those in the laboratory requires con-sideration of domains so wide as to exclude the recourse to fully resolved simulations.Using Gibson’s
C++ code
ChannelFlow , we scrutinize the effects of a controlledlowering of the numerical resolution on the decay of turbulence in plane Couetteflow at a quantitative level. We show that the number of Chebyshev polynomialsdescribing the cross-stream dependence can be drastically decreased while preserv-ing all the qualitative features of the solution. In particular, the oblique turbulentband regime experimentally observed in the upper part of the transitional range isextremely robust. In terms of Reynolds numbers, the resolution lowering is seen toyield a regular downward shift of the upper and lower thresholds R t and R g wherethe bands appear and break down. The study is illustrated with the results of twopreliminary experiments.Keywords: Plane Couette flow, Turbulence transition, Numerical simulationPACS: 47.20.FT 47.27.Cn The ‘laminar–turbulent’ transition in globally subcritical flows is far from being fullyunderstood. This is due to its abrupt and hysteretic nature, and to the fact that phasespace coexistence, typical of a subcritical bifurcation, has a nontrivial counterpart inphysical space, with laminar flow and turbulence coexisting in separate regions of theflow domain. Here we focus on plane Couette flow (PCF), the flow of a viscous fluid withkinematic viscosity ν sheared between two parallel plates at a distance 2 h , translating inopposite directions at speeds ± U . This flow configuration is free from global advection.The laminar flow is known to be linearly stable for all values of the Reynolds number R = U h/ν , whereas under usual conditions turbulent flow takes place for R large enough,typically R ∼ → turbulent’ (di-rect) and ‘turbulent → laminar’ (reverse). Many early studies have dealt with the directtransition, and especially with the dynamics of turbulent spots , by means of laboratoryexperiments or numerical simulations. More recently, experiments performed at Saclay [1]have shown that the reverse transition is marked by the occurrence of oblique turbulentbands , only observable in very large aspect ratio systems, in some range R g < R < R t .In the lowest part of this range, near R g ≃ R isfurther decreased below R g . Hysteresis is observed and sustained turbulent spots can beobtained by triggering the laminar flow with sufficiently large local perturbations when R > R g , whereas the laminar profile can be maintained up to much higher values of R provided that the experiment is sufficiently clean. At the upper end of the transitionalrange, the pattern disappears progressively and the transition from the turbulent bandsto featureless turbulence at R t ≃
410 is continuous. The term ‘featureless’ used here isborrowed from [2] where it served to describe the high- R turbulent regime beyond spi-ral turbulence in Taylor–Couette flow which corresponds to the oblique turbulent bandpattern in PCF.Besides laboratory experiments, numerical simulations of the Navier–Stokes equations(NSE) have provided invaluable information. An important output of early computationswas the concept of minimal flow unit (MFU) of size just necessary to maintain turbulencein a wall-bounded flow [3], a fundamental ingredient in the elucidation of the mechanismssustaining turbulence [4]. Later the MFU context was extensively used to study the decayof turbulence within the framework of dynamical systems theory [5]. Simultaneously,numerical simulations were also performed in wider domains, which lead to the discoveryof a large scale streamwise structures in turbulent PCF [6] and other wall-bounded flowsat Reynolds number somewhat beyond the transitional range defined above [7].Numerical studies specially dedicated to the problem of oblique turbulent bands arerecent. Soon after the experiments that put them at the forefront, Barkley & Tuckerman[8] succeeded in reproducing the fact by simulating the NSE in domains elongated in theexpected direction of the pattern’s wavevector but narrow in the complementary in-planedirection. These simulations gave useful information on the pattern, properly accountingfor the essential features of the laminar-turbulent alternation. The mechanism producingthe bands has however remained elusive up to now, and it is not clear whether periodicboundary conditions a few MFUs apart along the short dimension of these domains donot handicap our understanding of it. Although the occurrence of bands appears to bean extremely robust phenomenon, as our study will confirm, it thus seemed interestingto consider cases where the long-range streamwise coherence of the large scale streakystructures commonly observed in wall-bounded flows [7] was sufficiently well embraced.The coherence length of these structures being indeed at least one order of magnitudelarger than the streamwise length of the MFU, this revives simulations in large aspectratio, conventionally-oriented, domains. Such simulations again showed the occurrence ofoblique turbulent bands [9, 10].The computationally demanding character of these fully resolved numerical experi- The aspect ratio is the dimensionless size of the set-up, i.e. its lateral size in units of the cross-streamhalf-gap h . Subscript ‘g’ used hereafter stands for ‘global’ in the sense of ‘global stability threshold’, the value of R below which the laminar base flow profile is unconditionally recovered in the long term, i.e., whateverthe strength of the perturbation brought to the flow. y in termsof well-chosen ad hoc polynomials. The main characteristics of the transition were recov-ered at much lower numerical cost from a truncation of the expansion at first significantorder, which permitted simulations in very large aspect-ratio domains [12]. However, thetransitional range was lowered by a factor of two with respect to the experiments as aresult of insufficient energy transfer and dissipation in the cross-stream direction. Fur-thermore the oblique bands in the upper part of the transitional range were not obtained,presumably another effect of the lowered cross-stream resolution.The purpose of the work presented here is not to improve the model mentioned aboveby truncating it at higher orders, which is possible but very cumbersome and opaque, butto test this resolution effect in the context of direct numerical simulations, thus consid-ering the deliberate decrease of the spatiotemporal resolution as a systematic modellingstrategy. Our motivation is basically that, since qualitative and quantitative comparisonsof solutions obtained at different resolutions are easy, the degree of approximation canbe evaluated with some confidence. Having tested the reliability of this procedure, wemay expect to obtain clues on the mechanisms of band formation and decay directlyfrom the NSE at reduced numerical cost, in much the same way as lowering the size ofthe computational domain down to the dimensions of the MFU has helped toward theunderstanding of the self-sustaining process [13]. Finally, if quantitatively reliable lowresolution simulations can be performed, studying the statistics of the upper transitionat R t as well as the lower transition at R g will be possible in larger domains, duringlonger periods of time (the so-called ‘thermodynamic limit’ involved in analogies withthermodynamic phase transitions [15]), which will go in the same direction as in [12] butwithout the limitations of the model used in that work. Encouragement to follow theprogram sketched above can also be found in the work of Willis & Kerswell who obtainedenlightening results on statistical issues related to the dynamics of slugs and puffs forthe pipe flow transitional problem by modelling the flow through a drastic reduction ofthe azimuthal resolution [14]. In contrast with them, we shall also probe the reliabilityof the approach at a quantitative level by comparing results obtained when progressivelydecreasing the resolution progressively.Gibson’s well-known program ChannelFlow [16] is used throughout the presentstudy, taking for granted that it has already been abundantly validated and has served inmany high resolution studies. A first experiment, described in §
2, is focused on moderatelydeveloped turbulence somewhat above R t . The second experiment, in §
3, is devoted toa study of how the bifurcation diagram is damaged by a reduction of the spatiotemporalresolution in the transitional range: in a domain of size sufficient to contain one wavelengthof the band pattern, R is decreased by small steps from the turbulent regime down tothe laminar state and the values of R t and R g are systematically recorded. From thesetwo experiments, we determine an ‘optimal’ resolution above which the physics of thephenomenon is preserved, at the expense of a tolerable shift of these two thresholds.Some preliminary results supporting our prescription are presented in the last sectionending with a general discussion. 3 Decreasing the numerical resolution: qualitativeeffect on the turbulent state
Fields and functions in our
C++ code are defined as in
ChannelFlow . As to timeintegration, a backward formula is taken, which treats the viscous term implicitly andthe nonlinear term explicitly. The time step is adjusted so as to keep the CFL numberbelow 0.4 while being maintained smaller than 0.06 h/U in all our simulations. In thewall-normal direction, the spatial resolution is a function of the number N y of Chebyshevpolynomials used. The in-plane resolution depends on the numbers ( N x , N z ) of collocationpoints used in the evaluation of the nonlinear terms. From the 3/2 rule applied to removealiasing, this corresponds to solutions evaluated in Fourier space using N ′ x,z = 2 N x,z / δ eff x,z = L x,z /N ′ x,z = 3 L x,z / N x,z . In thefollowing, the resolution will everywhere be specified using the triplet ( N x , N y , N z ). In therange of Reynolds numbers of interest here (less than 500), simulations with N y ∼ L x,z /N x,z less than 0 . L x = L z = 32, lengths given in units of the half-gap h , see note 1) and R = 450, a value at which turbulence is permanently maintained. Initial conditions aretaken to be random and the simulations rapidly reach a steady state regime, typicallywithin less than 150 time units. Crude diagnostics about the flow regime are obtainedfrom an integral measure of the distance to the linear base flow derived from the quantity V − R ( u + v + w ) d x d y d z built in ChannelFlow , hereafter called ∆. ( V = 2 L x L z isthe volume of the full three-dimensional domain.) In order to get more local information,we have also considered the space-time dependance of the velocity departure u, v, w awayfrom the laminar base profile u b = y ˆx , and especially of the streamwise component u inthe plane y = 0, u ≡ u ( x, y = 0 , z ; t ), since it is a good tracer of any departure awayfrom the base state. As another observable, we have considered the perturbation kineticenergy E t = ( u + v + w ), either at a point, E t ( x, y, z ; t ), or its average over the gap,¯ E t ( x, z ; t ) = R +1 − E t ( x, y, z ; t ) d y . Fourier spectra of u or ¯ E t have been extensively used.Figure 1 displays the variation of the time average of ∆ over the interval t ∈ [400 , N y (left) and N x , N z (right) for a system of size L x × L z = 32 ×
32 at R = 450. Theexperiment where N y is varied (left panel) is performed with N x = N z = 128. Squaresmark the time-averaged ∆ and up/down triangles indicate the standard deviation ofits fluctuations around the average. The experiment with variable N x , N z (right panel)assumes N y = 21. It reports two cases. The first one is with N x = N z varying from 24 to256, where the squares indicate the time-averaged ∆ and up/down triangles the standarddeviation as before. The second one is with N z = 3 N x and the time-averaged ∆ markedwith stars, the corresponding standard deviation is the same order of magnitude and notshown for the sake of readability.A clear convergence of the results is observed as the resolution is increased, markedby a plateau reached for N y ∼
15 and N x = N z ∼
96. We are rather interested in theopposite limit of decreased resolution. Whereas a gentle trend is observed for N y ≥ N y = 9 and 7 behave differently and a divergence in finite time is even observedfor N y <
7. In fact, for the moderate Reynolds numbers of interest here, N y = 11appears to be a bound below which the numerical solution is no longer physical. Figure 24 N y ∆ N x =N z =128
96 644832 12824 160 192 224 2561440.250.270.290.310.330.35 ∆ N y =21 N z Figure 1: Variation of time-averaged quantity ∆ with the spatial resolution of the pseudo-spectral scheme, as measured by the number of modes ( N x , N y , N z ). Left: variable number N y of Chebyshev polynomials. Right: variable number of in-plane collocation points N x , N z . With L x = L z = 32, R = 450.displays the streamwise component of the disturbance in the plane y = 0. For N y = 9the solution shows an anomalously fine in-plane structure, as if turbulent energy couldonly be dissipated by being transferred to structures with fast variations in x and z .By contrast, the solutions obtained for N y ≥
11 display coarser velocity fluctuations,the case N y = 13 being already qualitatively similar to those for N y = 15, 21, and 25.These patterns are reminiscent of the very large streamwise streaky turbulent structuresobtained in wall-bounded flows [6, 7].These structures can further be identified in Figure 3 which presents Fourier spectraof u ( t = 2000) for the different resolutions considered. In all the graphs, the envelopesof the projections of the spectra are displayed, that is to say Σ( k x ) = max k z S ( k x , k z )and Σ( k z ) = max k x S ( k x , k z ), where S ( k x , k z ) = | ˆ u | is the Fourier spectrum of u . Theinformation contained in these quantities, though somehow limited, is more readable thanthe display of the full spectra while giving a proper account of the anisotropy of the flowthat clearly distinguishes the spanwise and streamwise directions.For variable N y , Fig. 3 (top), the cases N y = 7 and 9 are clearly not properly resolvedsince the envelopes of the projected spectra do not decay as k x and k z increase. Alreadypresent for N y = 11, the correct tendency strengthens as N y increases, though less rapidlyin the streamwise direction than in the spanwise direction. Along the k z axis, except for N y = 11, a clear peak is observable for k z ∼
7, i.e. λ z = L z / ∼ .
6; this value issomewhat larger than but of the same order of magnitude as the width ℓ z of the MFUmentioned above. In contrast, no peak at finite k x is observed in the streamwise directionbut a monotonic decrease as k x varies from zero to k x, max = N x /
3. These two featuresare in line with the observation of the very elongated streamwise streaky structures inwall-bounded flows. Similar trends are observed for N y = 21 and variable N x in Figure 3(bottom) where the spectra appear shifted with respect to one another due to the absenceof normalisation by the number of modes in order to improve the readability of the figure.The peak at k z ∼ N x,z ≥
64 but not at lower resolution which is not The units for k x,z are 2 π/L x,z . The maximum value of k x,z is N x,z / N ′ x,z = 2 N x,z / − N ′ x,z / N ′ x,z / x component of the velocity perturbationin the mid-plane y = 0 of full three dimensional snapshots of the numerical solutionsobtained in a domain 32 × ×
32 for R = 450 at t = 2000. N x = N z = 128 in-planecollocation points, and N y Chebyshev polynomials in the cross-stream direction, with N y ranging from 9 to 25. The x and z axes are respectively horizontal and vertical, like inall most of the other figures displaying patterns, except stated otherwise.6 k x
15 2519119713 k z k x k z Figure 3: Envelopes of Fourier spectra of the streamwise velocity component in the mid-plane u for L x = L z = 32, R = 450, and a single snapshot at t = 2000. Top: N y variable, N x = N z = 128. Bottom: N x = N z variable, N y = 21. Spectra are not scaled by the totalnumber of modes N x N z . For the top line N x N z is constant and spectra envelopes appearon top of each other, whereas for the bottom line the number of modes depend on theresolution so that they are shifted upward as it increases.7 −6 −5 −4 −3 −2 k x u −7 −6 −5 −4 −3 −2 k z u −7 −6 −5 −4 k x v −7 −6 −5 −4 k z v Figure 4: Fourier spectra of the streamwise ( u ) and wall-normal ( v ) corrections to thelaminar profile for the turbulent flow at R = 450, evaluated at y = 0 ( L x = L z =32, single snapshot at t = 2000, N y = 21 and N z = 3 N x variable). Top: streamwiseperturbation u . Bottom: wall-normal perturbation v . Left: streamwise component k x .Right: spanwise component k z . Successively: ( N x , N z ) = (24 , , , , , N x N z in order to show how well they pileup.surprising when we compare the corresponding wavelength L z / ≈ . δ eff x,z : for N x,z = 48, this makes δ eff x,z = 3 L x,z / N x,z = 1, hence less than fivepoints per spanwise wavelength, which really does not seem to be enough.Up to now we have not yet taken into account the fact that the velocity fluctuationsvary much more rapidly in the spanwise direction than in the streamwise direction. Herewe thus consider effective space steps systematically smaller by a factor of three along z than along x , i.e., N x = N z /
3. Results for ∆ are shown in the right panel of Figure 1 asstars. Except for the lowest value ( N x , N z ) = (24 , N x = N z when plotted as functions of N z , which means that in-planeresolution can be appreciated from the value of N z alone with N x down to a factor ofthree smaller. This observation is confirmed by the examination of the envelope spectraas defined above shown in figure 4. Both u and v ≡ v ( x, y = 0 , z ; t ) at t = 2000 havebeen considered. Here the Fourier amplitudes have been normalised by the correspondingnumbers of modes and it is clearly seen that they pile up for N x ≥ N z ≥ δ eff x = 3 L x / N x = 1 and δ eff z = 3 L z / N z = 0 . y < u + v + w > ( x , z ) −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 10 0.050.100.150.20 < u + v + w > ( x , z ) y
13 152125
Figure 5: Profiles of ¯∆( y ) as functions of the wall-normal resolution N y as interpolatedfrom values at collocation points using cubic splines. Left: lowest resolution, N y from 7to 13. Right: highest resolutions, N y from 13 to 25. ( L x = L z = 32, N x = N z = 128, R = 450, t = 2000.)How well the turbulent regime is reproduced in low resolution simulations has alsobeen tested by comparing in-plane averaged profiles of the streamwise velocity correction¯ u ( y ) = (1 /L x L z ) R u ( x, y, z ) d x d z and the squared distance to the base solution ¯∆( y ) =(1 /L x L z ) R ( u + v + w ) d x d z . Figure 5 displays the latter quantity for a series ofexperiments with N y varying from 7 to 25. Results for ¯ u ( y ) are similar. No averagingover time has been performed: a single snapshot at t = 2000 has been used in each case,which explains possible irregularities and a slight lack of y -symmetry for N y = 15 and21. From the figure, one can immediately see that profiles for N y = 7 and 9 are off, thatthose for N y = 11 and 13 progressively evolve so as to tend toward a limiting profile, andthat consistent results are obtained for N y ≥
15 in agreement with observations alreadymade when considering Figure 1 (left).
This preliminary study suggests that the featureless turbulent regime is reasonably wellrendered provided that N y ≥
15 and effective in-plane resolution δ eff x = 3 L x / N x ≤ δ eff z = 3 L z / N z ≤ .
33. It remains to test the effect of a lowering of the resolution on thepeculiarities of the ‘turbulent → laminar’ transition via oblique turbulent bands, which isdone in the present section.According to the experiments [1], the pattern can be characterised by a wavevector k patt with components k patt x = 2 π/λ x and k patt z = 2 π/λ z . Furthermore, λ x ≃
110 isobtained all along the transitional range and λ z varies from 50 for R = 395 close to thetop at R t ≃ R = 340 when the bands break down into patches that finallydecay below R g ≃ λ x and one spanwisewavelength λ z of the band pattern. Accordingly, we consider rectangular domains of size9
00 320 340 360 380 400 420 4400.160.180.200.220.240.260.280.300.32 ∆ R sustained t r an s i en t ≈ R g ≈ R t
11 13 15 21 33200250300350400450 laminaroblique bandsfeatureless
R N yR g R t Figure 6: Left: Time-averaged quantity ∆ as a function of R for N y = 21 with ( L x , L z ) =(110 ,
84) and ( N x , N z ) = (440 , R g and R t as functions of thenumber N y of Chebyshev polynomials for different in-plane resolutions and domain sizes,see text. L x × L z with L x ≥
110 and L z ≥
50. Except for N y = 7 and 9, oblique bands are obtainedwithout any difficulty in a full range of Reynolds numbers.Experiments are all done under the same protocol: a featureless turbulent regime isprepared at R = 450 and R is next decreased by steps. At every value of R , the simulationis pursued sufficiently long for the establishment of statistical equilibrium before furtherdecreasing R . Usually, this takes at least 5000 time units. The distance ∆ is recorded asa function of time. After elimination of a transient just after the change in R , its timeaverage and the rms value of its fluctuations are computed.The result of the experiment for N y = 21, ( L x , L z ) = (110 ,
84) and ( N x , N z ) =(440 , R , thetime averages of ∆ as squares and the standard deviations of fluctuations as up/downtriangles. (Estimates for transient states observed below R g are shown as open circles,see below.)A break in ∆ as a function of R is identified as the upper threshold R t . Visually,this break corresponds to the appearance of regions where the turbulence intensity isdepleted. A continuous diagonal band forms as R is decreased below R t , which, onceperiodically continued in the spanwise and streamwise directions features the expectedoblique pattern. The decrease of ∆ which measures the average distance to the laminarprofile is easily understood as the result of a decrease of the turbulent fraction, the ratioof the surface of the region still turbulent to the total surface of the system. This fractionis determined by thresholding the local mean perturbation kinetic energy ¯ E t ( x, z ; t ). Intheir simulations, Barkley & Tuckerman [8] observe that the turbulent intensity is largerinside turbulent bands than in the featureless regime. Quantity ∆ being an average overthe whole domain, by compensation this could lead to an underestimation of R t . Ratherthan trying to correct for this intensification effect we choose to define R t as given bythe position of the break further confirmed by visual inspection of the flow pattern.The second limit R g is determined as the value of R below which turbulence decays in Incidentally, turbulence intensification warrants study in domains of arbitrary shapes compatiblewith the band pattern. This phenomenon might indeed be specific to the narrow oblique domains usedin [8], since the periodic boundary forcing at short distance interferes with the large-scale streamwisecoherence of the streaks. R is observed between sustained banded turbulencewith limited fluctuations of ∆ and a turbulent regime bound to decay which displayslarge excursions toward values of ∆ well below the average. Below R g , taking the mean of∆ during the plateau that precedes the final decay produces the measurements displayedas open circles which appear to be well aligned with the other points corresponding tosustained turbulence.At any rate our aim here is not to perform a detailed statistical study of transientturbulence as in previous studies [5] but to locate R t and R g approximately as a functionof the resolution. A very conservative estimate of the precision with which R t and R g aredetermined is ±
5, which meets our purpose.The results of the most systematic study with ( L x , L z ) = (128 , N x , N z ) =(512 , R t , down-triangles for R g . Values at N y = 33, marked with hexagrams, are taken from the literature [9], in closeagreement with experimental observations [1]. The bifurcation diagram shows a regularshift of the interval in R where bands are observed. These results do not depend on theprecise size of the domain provided that it is sufficiently large to accommodate a band,as seen from a comparison of the results for L x = 128 and L z = 64 with those from theexperiment with L x = 110 and L z = 84 described above and reported as open circles at N y = 21. They are also not so sensitive to the in-plane resolution, as seen from the resultsat N y = 15 of the ‘high’-resolution experiment with ( N x , N z ) = 4 × ( L x , L z ) comparedto those of a ‘low’-resolution control experiment with N x = L x and N z = 3 L z displayedas open squares. A similar shift was observed by Willis and Kerswell in the pipe flowcase though in the opposite direction of increased thresholds [14]. However, in their case,the resolution decrease was in the azimuthal direction, which corresponds to our spanwisedirection. This difference is likely to play a role on the turbulence sustenance mechanisms,so that no definite inference can be made from this observation.Corroborating the observation of an anomalous behaviour of the fully turbulent statefor N y = 7 and 9 (Fig. 1, left; Fig. 2, top-left), banded patterns are not obtained in thesecases. On the contrary, a continuous decrease of ∆ is observed as R is decreased. For N y = 9, a featureless turbulent state similar to that at R = 450 (Figure 2, top-left) canindeed be maintained as low as R = 185 while for smaller R , a small-scale streamwisestructure happens to grow on top of a large scale spanwise modulation of the turbulenceintensity. This manifestly non-physical, slowly time-dependent, numerical solution finallydecays abruptly as the Reynolds number is further decreased but a detailed study of thistransition is pointless. Similarly, an even more exotic but equally uninteresting chaoticstate is obtained for N y = 7, and maintained at even lower R .A caveat is however required: at moderate (i.e. not the lowest) resolution, spuriousstable nontrivial nearly-steady states can be found at values of R in the lower part of therange where the bands exist and below, in the form of localised states resembling edgestates, e.g. [17]. An example is shown in Figure 7 (left). As understood from the displayof the local mean perturbation energy and the streamwise component of the perturbationvelocity, this solution has broken the symmetries of the plane Couette flow since it ispredominantly localised in y ∈ [0 ,
1] and not symmetrically distributed over [ − , R is varied. It turns out that this state exists in a limited range between a lower11ocal mean perturbation energystreamwise perturbation velocitycomponent u
12 13 14 15 16 17 18 19 20240250260270280290300310320 R sn N y Figure 7: Left: Spurious ‘edge state’-like numerical solution for N y = 19 close to thesaddle-node bifurcation point R ≈ . L x = 110, L z = 32, N x = 440, N z = 128); thissolution is mostly concentrated in the upper half of the channel ( y >
0) as understoodfrom the comparison of the variation ranges of vx ≡ u ( x, y, z, t ) for y = +0 .
766 and y = − .
766 at steady state ( t → ∞ ); Right: Variation of the saddle-node threshold R sn as a function of N y for the family of spurious numerical solutions.12addle-node threshold where it disappears and an upper threshold where it experiencesa Hopf bifurcation, and next nucleates a turbulent spot, which is expected and thus notstudied in detail. In fact, this solution belongs to a family that exists for a range of cross-stream resolutions N y . It was obtained from another one obtained during the decay of aturbulent state for N y = 15 and R = 270, then carried to N y = 13 on the one hand, and to N y = 17 and next N y = 19 on the other hand, but could be stabilised neither for N y = 11nor for N y = 21. Figure 7 (right) displays the threshold of the saddle-node bifurcationthrough which this solution family disappears. This threshold is seen to increase rapidlyand to lie below R g for N y = 13 or 15, and above for N y = 19, which is probably thereason why we cannot find it for N y = 21 without the help of a sophisticated continuationprocedure. This family of spurious solutions is most probably not unique. Unlike genuineedge states that are unstable by construction and proposed to represent gates on theboundary of attraction basin of the base flow, they are stable at least within some rangeof Reynolds numbers; the mechanisms by which they are maintained should however beessentially identical. We stress that this unwanted side effect of under-resolution here actson nontrivial but non-chaotic states and is not expected to spoil our study of turbulentregimes which are sensitive to noise inherent in chaos but statistically robust enough towithstand perturbations due to truncation errors.To conclude this section, in spite of the caveat above, evidence has been given that aslong as unsteady (turbulent) states are considered, reducing the wall-normal resolutionis a viable modelling strategy. Of course the best possible resolution is desirable but, inview of a quantitative reproduction of the transitional range of plane Couette flow, we canrecommend that the number N y of Chebyshev polynomials be kept larger than or equalto 11, whereas a tolerable shift of the upper and lower bounds of that range is obtainedprovided that N y ≥
15. In-plane resolutions with N x ≥ L x and N z ≥ L z also appearsatisfactory, with corresponding numbers of Fourier modes N x,z . We begin by presenting preliminary results obtained in parallel with the study presentedabove to illustrate how it can be used, before making more general comments on thenumerical approach to transitional wall-bounded flows. • Extreme case.
The emergence of turbulent bands first appears to be an extremelyrobust phenomenon. Here we show simulation results in the most extreme conditionsthat we have considered, namely N y = 11 and N x = L x , N z = L z . The lateral size ofthe domain is L x = 682 and L z = 381, which is already larger than in the early Saclayexperiments [18] and of the same order of magnitude as in the latest experiments [1] orthe simulations reported in [9] but with the computation power of a desk-top computer. From results in Fig. 6 (right), we expect R t ≃
270 and R g ≃ R = 310 ≫ R t . At time t = 500, R is suddenly reduced to Though we were unable to find any report on the existence of similar spurious numerical solutions—fragile with respect to a resolution change while fitting the framework of conventional bifurcation theory—in the open literature, we heard from Y. Duguet that the problem also arose in the search for exactsolutions in MFU-sized systems by Schmiegel (1999) and Gibson et al. more recently. Linux operated Dell computer with Intel Core2 DUO CPU E8400 3.00 GHz, 4 Go DDR2 memory.
500 1000 1500 2000 2500 3000 3500 4000 45000 0.050.100.150.200.25 ∆ t Figure 8: Quench experiment for N y = 11, N x = L x = 682, N z = L z = 341. Top: Timeseries of quantity ∆ for the values of R indicated. Bottom, from left to right and top tobottom: Snapshots of pattern obtained in the different asymptotic regimes obtained at t = 4000 after the quench for R = 250 (continuous bands), R = 240 (interrupted bands), R = 237 (border state), and at t = 500 for R = 235 (decaying spots).14ome final value R f . Time series of the quantity ∆ are shown in the top panel. Snapshotsof the solution are taken at regular times. The experiment is stopped at t = 4000 afterthe quench or when the laminar state is recovered.After the quench, a fast drop of ∆ is observed with a pronounced undershoot for R = 260, 250, and 240. During the drop, the evolution amounts to a general, moreor less uniform, breakdown of the turbulent state that, when the minimum is reached,leaves a few turbulent patches of limited extent. The turbulent fraction then re-increasesand bands form. For R = 260 and 250 continuous bands are obtained but, in contrastwith experiments where a single wave vector k patt was selected, here two wavevectorssymmetrical with respect to the streamwise direction ( ± k z ) dominate the pattern asdisplayed in the top-left image. The quench at R f = 240 ends with a regime where bandsare fragmented in the top-right image. For smaller R f , the undershoot is less pronouncedand the recovery stage slower. The asymptotic regime reached at R = 237 will be called a border state since—within the quench protocol—it sits on the frontier separating decayingfrom sustained regimes. As such, arguably, it could have been called an ‘edge state’ butthis term has been already extensively used in the study of the transition in terms ofdynamical systems [5] which focuses of the temporal dynamics of localised structures thatare exceptional limit sets. In contrast, the state shown in the bottom-left image appears typical of the spatiotemporal aspects of the transition. For R f = 230, the system does notrecover and the turbulent patches decay immediately. Apart from an important downwardshift of R t and R g , all of the observed scenario is identical to results described by Duguet et al. [9] obtained in a much better resolved context, or to Prigent’s experiments [1].From the above one could conclude that R g ≈ R g ≈ N x,z = 4 L x,z for results in that figure vs. N x,z = L x,z here. A similar effect was alreadyapparent for N y = 15 and the two resolutions considered: R g ≃
262 was obtained with N x,z = 4 L x,z (down triangle) and R g ≃
275 with N x = L x , N z = 3 L z (square), which isslightly better than now ( N z = L z ). Another, more likely, explanation can also be foundin the fact that the threshold obtained in the quench experiment is only an upper boundto R g : The quench protocol always involves a recovery stage starting at the end of theundershoot. The corresponding state is an assembly of turbulent spots similar to what ispictured in figure 8 (bottom, right) during the early decay at R = 235. In this respect,the quench experiment is an initial value problem similar to that of Duguet et al. [9], orto any experiment in which spots are triggered. In contrast, the adiabatic decrease of R isthe sole protocol supposed to yield, by continuation, the global stability threshold definedas the lowest value of R above which sustained turbulence has a non-empty attractionbasin. Except in a well-conducted adiabatic experiment, this attraction basin may indeedbe hard to find until one reaches values of R where it has a sizeable breadth, which issomewhat beyond R g .The superposition of two sets of bands illustrated above instead of one, as expectedfrom the experiments, is also the result of a resolution which, though qualitatively sat-isfactory, is nevertheless too poor since ongoing simulations with N y = 15 show a singledominant wavevector in the range of Reynolds numbers corresponding to bands. Interpret-ing the result in terms of pattern formation and Ginzburg–Landau envelope equations [1],when the resolution is low ( N y = 11) the coefficient of the cubic term | A ∓ | A ± accountingfor the interaction between wavevectors + k x and − k x (explicitly defined in the next para-graph) is under-estimated compared to the coefficient of the term | A ± | A ± accounting for15elf-interaction, so that each orientation can develop to form a rhombic pattern; a ratiocompatible with experimental findings is restored by increasing the resolution, ending ina stripe pattern ( N y = 15). On the other hand, it should be stressed that, despite theconstraints brought by the periodic boundary conditions, the pattern’s wavelengths andtheir variation are correctly reproduced since we can measure λ x ≈
114 = L x / L x / L x / R = 250 and 240 and λ z = 68 (= L z /
5) and85 (= L z /
4) for R = 250 and 240, respectively. • A better resolved experiment.
The case considered above stays too close to theboundary of the basin where DNS faithfully accounts for transitional plane Couette flow.The real motivation of our work is to consider less extreme simulation conditions in viewof a quantitatively better rendering of the different regimes observed. With our presentcomputational capabilities, this can only be done by increasing N y while simultaneouslykeeping a sufficiently high in-plane resolution at the expense of decreasing the size of thedomain. Accordingly we have considered a system of size sufficient to afford at least a fullelementary cell of the experimentally observed pattern, i.e. L x ∼ λ x and L z ∼ λ z . Theexample shown here is with L x = 110 and L z = 84 and a resolution N y = 15, N x = 512, N z = 256. It aims at the identification of an appropriate order parameter for the pattern[19]. The top row of Figure 9 displays snapshots of the local mean energy at differenttimes during a long simulation at R = 315. Band patterns with a single band leaning tothe right or to the left, with two bands leaning to the left, or a messy situation, can beobserved. Fourier analysis then yields either one dominant mode (well formed patterns,Fig. 9, left and centre) or several interacting modes (defective pattern, Fig. 9, right).In the present context, Fourier amplitudes are good candidates to play the role of orderparameters. They were indeed considered by Barkley et al. in [8] and [20]. The timeseries of the Fourier amplitudes shown as functions of time in Figure 9 (bottom) indicatean apparently random alternation of states that can be characterised by a dominantmode ( ± indicate left or right, and 1 or 2 the number of bands, i.e. k x = ± πn x /L x and k z = 2 πn z /L z , with n x = 1 and n z = 1 , R for different aspect-ratios L x,z . • Final remarks.
When combined with the reliability assessment provided in Sections 2and 3, these preliminary results bring an interesting contribution to the phenomenologyof plane Couette flow at minimal numerical cost. The pattern’s main characteristics arerecovered. Turbulent band formation appears to be a robust feature in the transitionalrange R ∈ [ R g , R t ]. The order of magnitude of the pattern’s wavelengths and theirvariations with R are also well reproduced in the simulations, even at the lowest possibleresolution. The price to be paid for the resolution decrease just seems to be a regularshift of the transitional range toward lower Reynolds numbers.These results also illustrate the specifically spatiotemporal features of the transitionalrange. In particular, turbulence decay may not well be rendered by the temporal approach16 −3 time (h/U) F ou r i e r a m p li t ude s ( A . U . ) (+1,1) (+1,2)(−1,2) (−1,2) (−1,2)(−1,1)(−1,1) (−1,1) Figure 9: Four snapshots of the solution for R = 315 ( N y = 15, L x = 110, N x = 512, L z = 84, N z = 256, x axis is vertical). Bottom: Time series of the amplitude of thedifferent dominant Fourier modes. 17mplemented in low-dimensional dynamical systems theory (chaotic transients [5]). Thelatter approach remains adapted to chaotic but spatially coherent dynamics in domainsa few MFUs wide. It can be of use to understand the nucleation of laminar patchesof limited extent within featureless turbulence. It is however unable to account for theregular regression of turbulent domains coexisting with laminar flow which marks thedecay stage in large aspect ratio systems, either in the laboratory or in the computer.This is attested by the variation of ∆ shown here for R = 235 and N y = 11 in Fig. 8, buttypical of better resolved cases.Interesting results have already be obtained in an elongated inclined domain [8, 20]but confinement by periodic boundary conditions in the direction parallel to the shortside perturbs the long-range streamwise coherence of the streaks (Fig. 2). This warrantsfurther study and leads us to suggest that one should consider domains that are at leastas large as one elementary cell ( λ x , λ z ) of the band pattern. As a modelling strategy within the extended systems perspective, results presentedhere drastically improves over the model previously elaborated in [11] and used in [12].Though that model is amenable to analytic treatment in view of the elucidation of themechanisms producing the experimentally observed turbulence modulation, the presentnumerical findings point out the role of its effective cross-stream resolution which is muchtoo low. On the other hand, the fact that limited resolution reproduces the main featuresof the dynamics of the flow in the transitional range suggests that the observed patternformation does not involve processes taking place in a thin boundary layer close to theplates where high-order cross-stream modes are of importance, but larger scale interactionsin the bulk of the shear where structures controlled by moderate-order modes operate (see[21] for the companion problem of puff sustainment in pipe flow). This observation couldmotivate the search for a model of intermediate complexity along the lines traced in [11]by pushing the Galerkin expansion a little further in view of an analytical approach.To conclude, plane Couette flow presents itself as an academic prototype of wall-bounded flows, with the extreme condition that it is linearly stable for all R , forcingits subcritical character. Low resolution simulations allows numerical experiments atminimal cost in circumstances where dynamics in physical space becomes more relevant(laminar/turbulent coexistence) than in phase space where the collection of competingexact solutions becomes increasingly large and their properties difficult to exploit. Ourunorthodox approach of decreasing the resolution in a controlled way may be consideredas a modelling methodology which will allow refined statistics, in view of Pomeau’s ther-modynamic analogy [15, 22], and help to pose appropriate questions about the physicsbehind band formation. The approach can certainly not be extended to the fully devel-oped regime nor to initial spot development where high resolution is an important issue,but moderately turbulent regimes involved in the transitional range of less academic flowscould take advantage of it when the lateral extension is of primordial interest.Acknowledgements: P.M. wants to thank Y. Duguet and G. Kawahara for interestingdiscussions related to this work, and the latter for his invitation to present it in Osaka.Helpful comments of a referee are also deeply acknowledged. This option was taken by Barkley in his latest work on band formation presented at the conference
New Trends On Growth And Form , Agay (France) June 20-25, 2010. eferences [1] Prigent A. et al. : Long-wavelength modulation of turbulent shear flows. Physica D , 100–113 (2003).[2] Andereck, C.D., Liu, S.S., Swinney, H.L.: Flow regimes in a circular Couette flowsystem with independently rotating cylinders. J. Fluid Mech. (2007).[8] Barkley, D., Tuckerman, L.S.: (a) Computational study of turbulent laminar pat-terns in Couette flow. Phys. Rev. Lett. (2009) 025301 [R]; 039904 [E].[13] Waleffe, F.: On a self-sustaining process in shear flows. Phys. Fluids (1986) 3-11. 1916] J. Gibson. .[17] Duguet, Y., Schlatter, Ph., Henningson, D.S.: Localized edge states in plane Couetteflow Phys. Fluid (2009) 111701.[18] (a) Bottin, S., Dauchot, O., Daviaud, F., Manneville, P.: Discontinuous transitionto spatiotemporal intermittency in plane Couette flow. Europhys. Lett. , 143–155 (1998).[19] Rolland, J., Manneville, P.: Oblique turbulent bands in plane Couette flow: Fromvisual to quantitative data, 16th International Couette–Taylor Workshop, PrincetonUniversity, September 2009. Ji Hantao, Smits, L. (eds.) [20] Tuckerman, L.S., Barkley, D., Dauchot, O.: Instability of uniform turbulent planeCouette flow: spectra, probability distribution functions and K − Ω closure model.Schlatter, Ph., Henningson, D.S. (eds.): 7th IUTAM Symposium on laminar-turbulent transition, Stockholm, June 2009. Springer (2010).[21] Shimizu, M., Kida, S.: A driving mechanism of a turbulent puff in pipe flow. FluidDyn. Res.41