Retraction: “Direct numerical simulation of the Ekman layer: A step in Reynolds number, and cautious support for a log law with a shifted origin” [Phys. Fluids 20, 101507 (2008)]

A retraction of our findings in Ref. 1 is necessary. The numerical accuracy was not up to the required standards, allowing an error of about 0.3 in $U+$ in the buffer layer, and was insufficient to support the logarithmic law that was proposed. New results on two levels of finer grids are presented and found to be close, and to compare very well with the channel direct numerical simulation (DNS) results of Hoyas and Jiménez.² Conversely, the excellent agreement with an analytical fit to experiments proposed by Chauhan et al.,³ well into the range of several hundred for $z+$ is lost.

The impetus for a deeper scrutiny of our results came from comparisons with other DNS data, which were not very close, even in the buffer layer, normally an area in which DNS is very reliable. The disagreement was exemplified at $z+=30$⁠, where the Ekman-DNS and Chauhan profiles were close to $U+=13$⁠, but high-quality channel DNS was near 13.3.² 

The resolution used for the Ekman layer was coarser than that in the channel. It was based on the guidelines of 1980s from the grid-resolution study by Spalart⁴ and will be called the 1988 resolution. Figure 2 in that article showed that a grid coarser than the 1988 resolution gave an underprediction by about two units of $U+$⁠, whereas a finer grid gave results in a cluster with all other runs (some of which used other domain sizes), with a scatter of about 0.3 units. This was considered adequate at the time, and these guidelines were followed in subsequent boundary-layer DNS studies, including all our previous Ekman work. While the conclusions in the 1988 article stand, the error bar on the velocity profile is larger than was thought then. For the purposes of the 2008 article, namely, determining the Karman constant precisely, via the derivative $dU+/dz+$⁠, this accuracy is not sufficient.

We conducted a new series of simulations to address this issue. The 2009 resolution amounts to using for each Reynolds number the grid previously used at the next higher Reynolds number. Recall that the steps were by a factor of $2$ for the Ekman Reynolds number Re, and by a factor of about 1.75 for the boundary-layer thickness in wall units, $δ+$⁠. The collocation points are now spaced in the wall-parallel directions by $Δx+=12$ and $Δy+=4$⁠, rather than about 20 and 7. Recall that with dealiasing, the shortest resolved wavelength is three times these spacings. In the vertical direction (also dealiased), the tenth collocation point from the surface is now at $z+=6$ or less, rather than 10, and, as the maximum CFL number is kept constant, the dimensionless time step $Δt+$ is reduced (for $Re=2000$⁠, from 0.165 to 0.129). The penalty for this improved resolution (which represents a substantial refinement for a spectral method) is that the value of the Ekman Reynolds number $Re=2828$ is now out of reach.

Tests were carried out to ascertain whether refinement in just one spatial direction, or in time resolution, would be sufficient. These were all unsuccessful, in that the results of the all-around fine grid were not reproduced. Refinement of two and three resolutions simultaneously were not tested, and so it remains possible that lower resolution in one or maybe two of $x$⁠, $y$⁠, $z$⁠, and $t$ might not adversely affect the quality of the simulation. The safer policy of refining by a similar factor in all directions was adopted.

Figure 1 compares results obtained using the 1988 and 2009 resolutions, and with a still finer grid (the previous resolution of $Re=2000$ applied at $Re=1000$⁠). As in Ref. 1, all curves are the result of averaging over horizontal planes and integer multiples of the inertial period $2π/f$⁠. Using these profiles to integrate the mean momentum equation recovers the mean surface shear stress to an accuracy of better than 1% in all cases.

At $z+=30$⁠, the profiles indeed rise from about 13 to about 13.3, thus matching channel DNS. This however seriously degrades the agreement with the experimental fit.³ The difference between the two finer grids is about 0.034, which then is our estimate of the residual uncertainty. The finest-grid profile is, again, slightly higher than the 2009 one. Figure 2 shows the Karman measure for the same simulations. Here, the difference between the two finest grids is about 0.006, compared to 0.027 between 1988 and 2009. This latter number is clearly excessive when attempting to select the Karman constant in the range from 0.38 to 0.42. The smaller number, 0.006, is probably adequate, but other obstacles remain for that exercise, as discussed below.

The “faintly circular” character of the previous study, mentioned in Ref. 1, was in fact more than faint: The study appeared to confirm the law of the wall, but only because the same wall-unit resolution was used throughout. All cases gave essentially the same profile, with essentially the same error. The other faintly circular aspect, related to the domain size, has not been tested.

Figure 3 is used to examine the agreement with channel DNS and the revised indications regarding the logarithmic law. The velocity profiles, not shown, make the 2009 curves indistinguishable from the channel curve up to roughly $z+=300$⁠, for Ekman $Re=2000$⁠. The Karman measure is more discriminating, yet indicates agreement within about 0.002 until the curves either give a strong appearance of departure from the law of the wall due to $z$ becoming commensurate with $δ$ (in other words, exiting the wall region) or begin to give an appearance of incomplete statistical convergence, revealed by a lack of smoothness. Our purely visual estimate of this location for the channel profile at $Reτ=2003$ is near $z+=250$⁠. At Ekman $Re=2000$⁠, the curve based on the velocity magnitude $Q$ is shown in addition to that based on the velocity $U$ in the direction of wall shear. It gives a further indication of nascent outer-layer (skewing) effects, creating a difference of 0.005 already at $z+=200$⁠, and somewhat worsens the failure to create a plateau, which remains the primary message of the figure.

The conclusion is that DNS of different flows agrees very closely where it can be expected to, except maybe for the channel flow at $Reτ=934$ (not shown). In that case, the Karman measure deviates from the $Reτ=2003$ result by 0.01 near $z+=105$⁠, and later exceeds 0.5 (thus clearly entering the outer layer) around $z+=750$⁠. The ratio between the two, namely, 0.14, is much smaller than for the Ekman layer at $Re=1000$⁠, in which the numbers are about 190 and 350 and the ratio 0.54. Thus, the statement in Ref. 1 that “for reasons unknown,” channel flow does not establish the law of the wall as fast as the Ekman layer, remains correct. It is possible that the Ekman layer experiences a fortuitous cancellation between upward or downward effects when the outer-layer dynamics begin to alter the velocity profile.

Compared to Fig. 6a in Ref. 1, the impact of the finer grid on the Karman measure is strong. Instead of peaking at about 0.411, it now peaks at 0.434 and ebbs below 0.385, for $Re=2000$ (0.375 if using $Q$⁠). Lacking the $Re=2828$ case, it is more difficult to estimate where the outer layer begins. In addition, proposing a shifted origin as we did before is now less convincing because of the much steeper slope the Karman measure displays in the [100,200] range; to compensate for it, the $a+$ shift would have to be much larger than 7.5, and therefore more difficult to reconcile with the thickness of the buffer layer. The much larger distance between the DNS and the experimental fit also damages the confidence that grew from the agreement of two such independent sources. This remains to be resolved.

Finally, the similarity theory described in Ref. 1 is affected in two ways: First, it assumes that a meaningful logarithmic law exists at the Reynolds number at which it is applied, which is evidently not the case for the range covered by the new results (again, we are reluctant to accept the large origin shift which would now be needed). Second, the limited extent of this range means that it is not at present possible to determine the constants with satisfactory precision simply by using the three available data points for the variables $u∗$ (surface friction velocity) and $α0$ (surface shear angle).

The unfortunate conclusion from this new study and Ref. 2 is that the goal of determining the constants of the logarithmic law from DNS, which is fully dependent on finding a plateau in Fig. 3, has receded by years. It does not appear safe to predict that this will be possible with less than a factor of about 5 in Reynolds number, and therefore a factor of near 600 in computing cost, or about 15 years at Moore’s rate. In the late 1980s, it was expected that the increase in Reynolds number $δ+$ by an order of magnitude that has now been achieved (with essentially the same codes) would ensure a clear message. Instead, DNS as well as experiments in pipes and boundary layers indicate that the entry into the log law is incomplete up to several hundreds in $z+$⁠, and even occurs from above instead of from below, as appeared natural. Other surprises are possible.

Some of this work was done as part of the UK Turbulence Consortium, sponsored by the Engineering and Physical Sciences Research Council (Grant No. EP/D044073/1). Computations were made on the U.K./EPSRC HPCx system.