Turbulent Couette flow up to ${{Re}}_\tau =2000$

Abstract Two simulations of turbulent Couette flows were performed at friction Reynolds numbers of 1000 and 2000 in a large box of dimensions $L_x=16{\rm \pi} h$, $L_y=2h$ and $L_z=6{\rm \pi} h$, where h is the semi-height of the channel. The study focuses on the differences in the intensity and scaling of turbulence at these two Reynolds numbers. The 2000 case showed a lack of a clear log layer with a higher value of the Von Kármán constant $\kappa$ than Poiseuille channels. The intensities were well-scaled in the buffer layer and below, with a second maximum of the streamwise intensity at approximately 350 wall units. Contrary to Poiseuille channels, the dissipation scales close to the wall in wall units. This fact can be attributed to the constant value of the derivative of the streamwise intensity in wall units. The intensities of the 2000 case showed remarkable differences compared with those at Reynolds number 1000 at the channel centre, likely due to the organization of large scales of the streamwise fluctuactions, $u$. These large scales were thought to be considered ‘infinite’. However, for the 2000 case, while all the structures have a width of $\ell _z \approx 6/8{\rm \pi} h$, their length varies from $\ell _x \approx 6{\rm \pi} h$ to $\ell _x \approx 16{\rm \pi} h$, which clearly contradicts the trends obtained in the past. This is a new effect that has not been reported for turbulent Couette flow and points to the uncertainty and sensitivity that is observed for certain statistical quantities.


Introduction
Direct numerical simulation (DNS) has become one of the main tools in studying wall turbulence.This fact is especially true in the case of Couette flows, which are investigated in the present study; see figure 1.The flow is confined between two parallel plates.The bottom plate remains stationary while the top one moves at a velocity U ∞ .As in pressure-driven channels, i.e.Poiseuille flows, the main control parameter is the friction Reynolds number Re τ = u τ h/ν, where u τ is the friction velocity, ν is the kinematic viscosity, and h is the channel half-width.Since the seminal work of Kim, Moin & Moser (1987), the turbulence community focuses on plane Poiseuille flows.The value of Re τ has increased from 180 to 10 000 in 35 years (Oberlack et al. 2022).However, considerably less attention has been paid to plane Couette flows.Several reasons make these flows harder to study.From an experimental point of view, the set-up with a moving wall is difficult to implement properly (Tillmark 1995).The main issue for numerical studies is the long and wide structures in the form of streamwise rolls that exist in turbulent Couette flows.Such structures are coherent regions of either positive or negative streamwise velocity fluctuations.They have been observed both experimentally (Tillmark 1995;Kitoh, Nakabyashi & Nishimura 2005) and numerically (Bech et al. 1995;Komminaho, Lundbladh & Johansson 1996;Tsukahara, Kawamura & Shingai 2006;Pirozzoli, Bernardini & Orlandi 2011;Bernardini, Pirozzoli & Orlandi 2013;Avsarkisov et al. 2014;Pirozzoli, Bernardini & Orlandi 2014), persisting even if the flow is influenced by a cross-flow (Kraheberger, Hoyas & Oberlack 2018) or a streamwise pressure gradient (Gandía-Barberá et al. 2018).Lee & Moser (2018) obtained structures at least as long as 310h for Re τ = 500.On the other hand, Gandía Barberá et al. (2021) has shown that stratification processes in thermoconvective flows can effectively remove them.
In a recently published theoretical work by Dokoza & Oberlack (2023) on the roll-cell structures in turbulent Couette flow, these structures are analysed for large longitudinal extension (α → 0), where α represents the wavenumber in streamwise direction, in the limit of large Reynolds numbers.Therein, based on the resolvent analysis, the authors show that in the distinguished limit α → 0 and Re → ∞, where Re here, unlike Re τ , is based on the bulk velocity, the product Reα = O(1) delivers invariant structures, and, as a result, the rolls grow with increasing Reynolds number.The original idea goes back to an analysis by Yalcin, Turkac & Oberlack (2021), who showed the corresponding invariance for the Orr-Sommerfeld equation.
For Poiseuille flow, Lozano-Durán & Jiménez (2014) noticed that even relatively small computational boxes of streamwise and spanwise sizes of only 2πh × πh can satisfactorily recover the one-point statistics of the flow.Lluesma-Rodríguez, Hoyas & Peréz-Quiles (2018) also proved a similar result, considering the energy equation.However, this is a controversial point in the Couette flow case.As mentioned, the long and wide structures can be as long as the computational box itself (Lee & Moser 2018).Consequently, a wide range of dimensions of the computational box has been studied.relevant publications of Couette flow studies, together with Re τ , the computational box dimensions and the mesh resolution.
The most crucial point about these rolls is that they can significantly influence the flow statistics (Alcántara-Ávila et al. 2019).In that work, the authors performed several simulations of a turbulent Couette flow at Re τ = 180 with different spanwise sizes, ranging from π/2h to 6πh.They found that the number of velocity rolls does not affect the mean flow.However, intensities are greatly affected, mainly in the centre of the channel.Thus, a wide computational box that can accommodate several rolls is needed to simulate the flow accurately.
Our present study aims to analyse the kinematics of Couette flow at a larger Reynolds number, up to Re τ = 2000 in a wide and long box.We found that it is possible that the rolls are not as large as previously thought.Roll length can be again finite for this friction Reynolds number.Apart from that, particular emphasis will be given to several points raised recently in Poiseuille flows, including the value of the von Kármán constant, the scaling of the intensity of the streamwise velocity and the scaling of the dissipation (Hoyas et al. 2022).

Numerical methods and simulations
This work presents two new simulations for friction Reynolds numbers Re τ = 1000 and Re τ = 2000.Both computations have been performed within a computational box of L x = 16πh, L y = 2h and L z = 6πh with spanwise and streamwise periodicities.The streamwise, wall-normal and spanwise coordinates are x, y and z and the corresponding instantaneous velocity components are U, V and W, respectively.Defining the space-time average operator • x i as (2.1) The value of φ x can be thought as the mean in x of the time-averaged field φ.Statistically averaged quantities in x and z are denoted by an overbar, φ = φ xz , whereas fluctuating quantities are denoted by lowercase letters, i.e.U = U xz + u = Ū + u.Primes are reserved for intensities: u = uu 1/2 .
The numerical methods of the two new DNS are based on those described by Kim et al. (1987).A similar algorithm, including the energy equation, is described in  Avsarkisov et al. (2014).This code is highly parallelizable, including the input/output operations, allowing high-resolution simulations (Yousif et al. 2023) or high-demanding computational power (Hoyas et al. 2022).
Table 2 summarizes the parameters of all the Couette flow simulations used in this paper.Information about the two new cases described above, namely C1000 and C2000, are given in the first two rows.The wall-normal grid spacing of these two cases is adjusted to keep the resolution Δy ∼ = 1.5η approximately constant in terms of the local Kolmogorov scale η = (ν 3 /ε) 1/4 .Using wall units, Δy + varies from 0.32 at the wall to Δy + 8.9 at the centreline.
The running times given in table 2 are provided regarding eddy turnover time units, u τ T/h.The initial file of these simulations was taken from a smaller Reynolds number simulation.In particular, the initial file for the C2000 simulation was taken from the C1000.The table does not include the transition until a statistically steady state is reached.One of the tools used to assert whether such a statistically steady state has been reached is to compute the total shear stress, figure 2(a).This stress should equal 1 in Couette flows, as integration and non-dimensionalization of the x-component of the mean momentum equation (Avsarkisov et al. 2014) yields The parameter ε is defined as the L2 norm of the difference between the two sides of (2.2) (Vinuesa et al. 2016).This difference is below 3 × 10 −3 in all cases, as shown in figure 2(b).The ε value in table 2 indicates that the simulation statistics have converged.Notice that every simulation follows this convergence criterion, due to the mesh sizes being similar.This is remarkable as Pirozzoli et al. (2014), and Lee & Moser (2018) use a different method from LISO.
We have also computed a bound for the error of the statistics as in Hoyas & Jiménez (2008).While the standard deviation of the data is below 0.4 % for Ū, it is considerably larger, up to 7 %, for the intensities.As explained later, this is an effect of the very large structures of the flow.As the number of samples is 85 000, the standard error is below 0.1 % in every statistic presented.

Mean velocity
The mean velocity profile is shown in figure 3 The thin solid line represents the classical logarithmic law of the wall for a von Kármán constant κ = 0.397, the value obtained for the C2000 simulation.For the case C1000, a value of κ = 0.41 fits the data a bit better, while in previous Couette simulations, the following values were obtained: κ = 0.41 (Pirozzoli et al. 2014), κ = 0.41 (Avsarkisov et al. 2014) and κ = 0.384 (Lee & Moser 2018).A constant region, such as that appearing in Hoyas et al. (2022), is not observed for any cases, including C2000.As expected (Lee & Moser 2015), even the highest Reynolds number presented here seems too low to produce This value has been discussed in several works.Reichardt (1959) and Busse (1970) determined a value of Ψ = 0.25 at infinite Reynolds number, while Lund & Bush (1980) performed an asymptotic analysis and concluded that it approaches zero as Re → ∞.In our case, however, we find Ψ 1000 = 0.066 and Ψ 2000 = 0.097, which are in line with the results of Lee & Moser (2018).Hence, as Ψ 2000 > Ψ 1000 , this may indicate a change of tendency or a hint of the difficulty of accurately computing statistics on the outer flow due to the existence of very large structures, see figure 3(b).This figure shows the value of Ψ evaluated after averaging 100 images of the flow.The mean of Ψ is basically the same as that computed with the whole database, but the differences across z are very high, probably due to the large structures of the flow.Finally, as Ψ involves wall-normal derivatives, a very long time for computing accurate statistics should be expected (Hoyas et al. 2024).

Velocity fluctuations
Neither v + or w + collapses exactly in wall units, see figure 4(a).As in Poiseuille flows (Hoyas et al. 2022), the maximum of these intensities is still growing with Re τ .The value of v + is remarkably similar to that predicted for Poiseuille flows in Spalart & Abe (2021).However, there seems to be a greater dispersion in the value in the channel centre figure 4(b), with the C2000 case being rather different.A simulation at a larger Reynolds number seems necessary to elucidate if this is a change of tendency or, again, an effect of the structures (Alcántara-Ávila, Hoyas & Pérez-Quiles 2018).
The intensity of u, u , is shown in figure 5(a).Close to the wall, the peak of the u + does not seem to increase appreciably with the Reynolds number but reaches a finite value (last column of table 2).These numbers are appreciably larger than that of Hoyas et al. (2022) for Re τ = 10 000, max(u + ) = 3.07.Several experimental (Hultmark et al. 2012;Willert et al. 2017) and theoretical studies (Chen & Sreenivasan 2021) suggested that this limit could be bounded in wall turbulence.Furthermore, this limit does not depend on the channel length, as Lee and Moser (2018) found for Re τ = 500 and L x = 100πh.This boundness is a remarkable difference between Couette and pressure-driven flows, as this maximum grows without bound in the range of Reynolds numbers described up to now (Marusic, Baars & Hutchins 2017;Hoyas et al. 2022).The second maximum in u + , at approximately y + = 270, appears at a Reynolds number far smaller than in Poiseuille channels (Hoyas et al. 2022) or boundary layers (Samie et al. 2018).
The agreement between Couette and Poiseuille flows is remarkably good in wide regions of the turbulent budget terms of u i u i : where D/Dt is the mean substantial derivative.The terms on the right-hand side are production, dissipation, turbulent diffusion, pressure strain, pressure diffusion and viscous diffusion.The exact definition of these terms can be found in Hoyas & Jiménez (2008) and Avsarkisov et al. (2014).For the sake of brevity, we will only discuss B 11 here.For the other cases, see the availability statement about the data.Figure 5(b) presents data for Couette (solid lines) and Poiseuille (dashed lines).The Couette cases follow the same colour code as in table 2. The two Poiseuille channels are at Re τ = 2000 (magenta) (Hoyas & Jiménez 2006) and Re τ = 10 000 (purple) (Hoyas et al. 2022).The first conclusion is that the scaling is perfect for all cases above y + ≈ 10.The peak of the production term at y + ≈ 12 can be considered universal.
A second point is the anomalous scaling of the dissipation (Hoyas & Jiménez 2008).As it is well known, the dissipation does not collapse in wall units for the Poiseuille cases (Hoyas et al. 2022).This situation changes for Couette flow.As shown in figure 5(b), the dissipation value at the wall is almost constant, with a value of (d/dy)u + | y + =0 = 0.4665 and 0.4653 for C1000 and C2000, respectively.This behaviour can be linked to the finite value of the maximum of u + (Alcántara-Ávila, Hoyas & Pérez-Quiles 2021;Hoyas et al. 2022).

Superstructures in Couette flows
As stated above, the most remarkable difference between Couette and Poiseuille flows is the existence of very long rolls with a width of approximately δ ≈ 6/8πh.
The footprints of these structures fill the whole channel, as shown in figure 6.Here, we have plotted the x-average of the streamwise fluctuation, u x , as the background colour.The value of the vector field ( v x , w x ) is represented with arrows.As we can see, it is possible to identify four structures with a width of z ≈ 6/8πh, approximately.This width is similar to the previous ones reported (Lee & Moser 2018).The length of these structures can be obtained by studying the correlation in x of u: This correlation is shown at the channel centre, figure 7(a,b).The white lines mark the boundary among the rolls.Lee & Moser (2018) showed that for Re τ = 500, the structures were extremely long, with a length of at least 100h.Alcántara-Ávila et al. (2019) obtained a similar result for thermal Couette flow.Our results offer a different trend than the previous results at smaller Reynolds numbers.In the C1000 case, the correlation is large for the four structures, filling the entire box length, while the structures of the C2000 case from z/h = 0 to 3π are smaller than the box.This can be appreciated better in figure 7(c).
Here, the value of R x can be seen at x = 6πh.The maximum of this correlation is always greater than 0.3 for the C1000 case, but it is negative in some z/h points for the C2000 case, indicating that the length of the coherent structure is less than x .Therefore, some parts of the channel can keep very large rolls, while in other regions, these rolls are limited to lengths below x = 6πh.This is a new effect that has not been reported for turbulent Couette flow and points to the uncertainty that is observed for certain statistical quantities, such as the slope parameter ψ in (3.1).

Conclusions
This paper presents two simulations at two different friction Reynolds numbers, namely 1000 and 2000.The simulations were carried out in a very large box of size L x = 16πh, L y = 2h and L z = 6πh.For the C2000 case, a clear log layer is only observed in parts.Through a fitting, a value of κ = 0.397 has been obtained.This value is a bit higher than Poiseuille channels and similar to other Couette ones.The scaling of the intensities is perfect in the buffer layer and below.The first maximum of u + does not grow anymore, and a second maximum has clearly appeared at y + ≈ 350.This value is similar to that expected for Poiseuille channels (Hoyas et al. 2022), but at a far lower Reynolds number.
The scaling failure of the dissipation is also not present anymore.This can be traced to the constant value of (d/dy)u + at the wall.The intensities of the C2000 case are remarkably different from those at lower Reynolds numbers at the channel centre.This is probably due to the organization of the large scales of u.These large scales were thought to be almost infinite.We obtained the same result for the C1000 case.However, for the C2000 case, some rolls present a finite size with a length of approximately x = 6π.The width of every roll is approximately z = 6/8π.Simulations at a larger Reynolds number are needed to elucidate if these structures are bounded.
Funding.The authors gratefully acknowledge the provision of computing time on the Gauss Centre for Supercomputing e.V. on the GCS Supercomputer SuperMUC-NG at Leibniz Supercomputing Centre under project number pr92la.M.O.acknowledges partial funding by the German Research Foundation (DFG) under the project OB 96/48-1.S.H. is partially funded by project PID2021-128676OB-I00 by Ministerio de Ciencia, innovación y Universidades/FEDER. Declaration of interests.The authors report no conflict of interest.

Figure 1 .
Figure 1.Sketch of plane Couette flow driven by the velocity U w of the upper wall.The dimensions of the computational box are L x × 2h × L z .

Figure 2 .Figure 3 .
Figure 2. Lines as in table 2. (a) Shear and Reynolds stress: d Ū+ /dy + (dash-dot), uv + (solid).Symbols as given in table 2. Notice how the two cases at Re τ = 1000 are basically the same.(b) Difference between left-hand side and right-hand side of (2.2).
(a)  in terms of the indicator function, Γ = y + ∂ y + Ū+ .This function shows a plateau if the scaling for the logarithmic layer Ū+ = κ −1 log( y + ) + B holds, where κ is the von Kármán constant(Oberlack et al. 2022) and B is the interception constant.

Figure 4 .
Figure 4. Lines as in table 2. Thin lines: v + .Thick lines: w + .The y scale is wall units y + (a) and outer units y/h (b).

Figure 6 .
Figure 6.The x-averaged field u x .The value of the vector field ( v x , w x ) is represented with arrows.(a) C1000 case.(b) C2000 case.

Figure 7 .
Figure 7. (a,b) R x (s) for C1000 (a) and C2000 (b).In both figures, we have added φ x , showing the size in z of the pair of rolls.(c) R x (s) at x = 6πh.The different rolls are marked using a dashed vertical line.

Table 1 .
Table1contains the most Most relevant publications of DNS of plane Couette flow.The second column shows the ranges of Re τ simulated for the dimensions of the computational box shown in columns three and four.The last three columns show the mesh resolution.

Table 2 .
Lee & Moser (2018)2014)ons discussed.Here L x and L z are the numerical box's periodic streamwise and spanwise dimensions; h is the channel half-height; x + and z + are the resolutions in terms of dealiased Fourier modes; N x , N y , N z are the numbers of collocation points.The time span of the simulation is given in terms of eddy turnovers u τ T/h, where T is the computational time; ε is a measure of convergence, defined inVinuesa et al. (2016).The maximum of u + is given in the last column.Line styles given in the second column are used throughout the paper.Data references are C986Pirozzoli et al. (2014),C550 Avsarkisov et al. (2014)and C500Lee & Moser (2018).