Structure and interpolation of the turbulent velocity profile in parallel flow

TL;DR

This paper examines the classical turbulence velocity profile theory, identifies discrepancies with experimental data, and proposes a correction based on pressure gradient effects to improve the theoretical description.

Contribution

It introduces an analytical correction to the classical turbulence profile theory accounting for pressure gradient effects, enhancing its accuracy.

Findings

Discrepancies in the logarithmic layer can be explained by pressure gradient effects.

A universal correction improves the match between theory and experiments.

The correction enables better analytical interpolation of velocity profiles.

Abstract

The classical scaling theory of turbulent parallel flow provides a framework for the description of the mean velocity profile through two functions of one variable, traditionally named law of the wall and law of the wake, and a universal logarithmic law characterized by von K\'arm\'an's constant. Despite its widespread adoption in research and in teaching, discrepancies between this theory and both experiments and numerical simulations have been repeatedly observed in the literature. Recently we have shown that in the logarithmic layer such discrepancies can be physically interpreted and analytically accounted for through an equally universal correction caused by the pressure gradient. This finding opens the way to a likewise improvement in the description of the law of the wall and of the law of the wake, an analytical interpolation of either of which is often useful for practical applications.

Tables

Table 1. Table 1: Curve labelling and data sources

Acronym	Geometry	${\mathrm{Re}}_{\tau }$	Technique	Reference
HJ550	plane	547	DNS	[9]
HJ1E3	plane	934	DNS	[9]
HJ2E3	plane	2004	DNS	[9]
DJ4E3	plane	4179	DNS	[11]
LM550	plane	543	DNS	[16]
LM1E3	plane	1001	DNS	[16]
LM2E3	plane	1995	DNS	[16]
LM5E3	plane	5186	DNS	[16]
BPO550	plane	550	DNS	[5]
BPO1E3	plane	999	DNS	[5]
BPO2E3	plane	2021	DNS	[5]
BPO4E3	plane	4079	DNS	[5]
Kea360	pipe	361	DNS	[10]
Kea550	pipe	550	DNS	[10]
Kea1E3	pipe	999	DNS	[10]
WM700	pipe	685	DNS	[30]
WM1E3	pipe	1142	DNS	[30]
N300	pipe	294	Exp.	[21]
N700	pipe	657	Exp.	[21]
N2500	pipe	2511	Exp.	[21]
N4500	pipe	4537	Exp.	[21]
N14000	pipe	14200	Exp.	[21]
N35000	pipe	35200	Exp.	[21]
N55000	pipe	55500	Exp.	[21]
Sea2E3	pipe	1985	Exp.	[26]
Sea3E3	pipe	3334	Exp.	[26]
Sea5E3	pipe	5411	Exp.	[26]
Sea1E4	pipe	10480	Exp.	[26]
Sea2E4	pipe	20250	Exp.	[26]
Sea4E4	pipe	37690	Exp.	[26]
Sea7E4	pipe	68160	Exp.	[26]
Sea1E5	pipe	98190	Exp.	[26]
PBO260	couette	260	DNS	[24]
PBO500	couette	507	DNS	[24]
PBO1E3	couette	986	DNS	[24]

Table 2. Table 2: Coefficient of the linear velocity correction as fitted to individual velocity profiles.

Reference	${\mathrm{Re}}_{\tau }$	$A_0$	$A_1$
[9]	547	1/0.392	1.01
[9]	934	1/0.392	1.10
[9]	2003	1/0.392	1.00
[11]	4179	1/0.392	0.76
[16]	543	1/0.392	1.05
[16]	1000	1/0.392	1.06
[16]	1995	1/0.392	1.13
[16]	5186	1/0.392	0.84

Table 3. Table 3: Best fit coefficient of the linear velocity correction for pipe flow.

Reference	${\mathrm{Re}}_{\tau }$	$A_0$	$g{A}_1$
[10]	361	1/0.392	1.83
[10]	550	1/0.392	1.91
[10]	999	1/0.392	2.04
[30]	685	1/0.392	2.22
[30]	1142	1/0.392	2.43

Table 4. Table 4: Coefficients of the exponential-Padé approximant of the law of the wall.

Equations

$$-p_x=4{\tau }_w/D_H$$

$$u^{+}=u^{+}(z^{+};h^{+})$$

$$u^{+}\simeq f(z^{+});$$

$$u^{+}\simeq u^{+}(h^{+})-F(Z);$$

$$u^{+}\simeq A\log (z^{+})+B=A\log (Z)+A\log (h^{+})+B,$$

$$u^{+}=f(z^{+})+u(h^{+})-F(Z)-A\log (Z)-A\log (h^{+})-B,$$

$$u^{+}=f(z^{+})+G(Z).$$

$$G(Z)=-F(Z)-A\log (Z)+C,$$

$$u^{+}(h^{+})={\kappa }^{-1}\log (h^{+})+B+C.$$

$${\kappa }^{-1}=[u^{+}(h^{+})-u^{+}(50)]/[\log (h^{+})-\log (50)];$$

$$\begin{gathered} u^{+}=A_0\log (z^{+})+A_{1}g{\mathrm{Re}}_{\tau }^{-1}{z}^{+}+B= \\ =A_0\log (Z)+A_{1}gZ+A_0\log ({\mathrm{Re}}_{\tau })+B, \end{gathered}$$

$$u^{+}=f_0(z^{+})+{h^{+}}^{-1}{f}_1(z^{+})+{h^{+}}^{-2}{f}_2(z^{+})+\dots ,$$

$$u^{+}(h^{+})-u^{+}=F_0(Z)+{h^{+}}^{-1}{F}_1(Z)+{h^{+}}^{-2}{F}_2(Z)+\dots .$$

$$\begin{gathered} u^{+}=f_0(z^{+})+{h^{+}}^{-1}{f}_1(z^{+})+{h^{+}}^{-2}{f}_2(z^{+})+\dots + \\ +G_0(Z)+{h^{+}}^{-1}{G}_1(Z)+{h^{+}}^{-2}{G}_2(Z)+\dots \end{gathered}$$

$$u^{+}=f_0(z^{+})+{h^{+}}^{-1}\left[A_{1}g{z}^{+}+f_1(z^{+})\right]+\dots +H_0(Z)+\dots$$

$$H_0(Z)=G_0(Z)-A_{1}gZ.$$

$$u_1^{+}(z^{+})=f(z^{+})+G(z^{+}/h_1^{+})$$

$$u_2^{+}(z^{+})=f(z^{+})+G(z^{+}/h_2^{+})$$

$$\Delta G_{12}(z^{+}/h_1^{+})=u_1^{+}(z^{+})-u_2^{+}(z^{+}),$$

$$\Delta G_{12}(Z)=G(Z)-G(mZ)$$

$$G(Z)=\sum _{n=0}^{\infty }\Delta G_{12}\left(m^{n}Z\right).$$

$$G(mZ)=\sum _{n=1}^{\infty }\Delta G_{12}\left(m^{n}Z\right),$$

$$f(z^{+})=u_2^{+}(z^{+})-G(z^{+}/h_2^{+}).$$

$$\Delta f_{21}(Z)=u_2^{+}(h_2^{+}Z)-u_1^{+}(h_1^{+}Z)=f(h_2^{+}Z)-f(h_1^{+}Z).$$

$$\Delta f_{21}(Z)\simeq A\log (h_2^{+}/h_1^{+})\text{for }1/h_1^{+}\ll Z\le 1;$$

$$u^{+}=\log (z^{+})/0.392+4.48,$$

$$H(Z)\simeq -0.57Z^7.$$

$$G(Z)\simeq A_{1}Z-0.57Z^7$$

$$H(Z)=-0.67Z^7.$$

$$G(Z)=g{A}_{1}Z-0.67Z^7$$

$$H(Z)=\frac{x-0.5}{\exp (-25(x-0.5))-1}$$

Full text

††journal: European Journal of Mechanics B/Fluids

1 Background and introduction

The classical scaling theory of the turbulent velocity profile, based after Prandtl, von Kármán and Millikan on the distinction between a wall layer and a defect layer, is a mainstay of turbulence theory. It underlies the friction law of pipes and channels and also constitutes a basic ingredient of the description of the turbulent boundary layer. Therefore, the extraction of the empirical coefficients (and functions, such as the universal law of the wall and Coles’ law of the wake) contained in this theory from experimental data and numerical simulations is a continuing activity that absorbs the attention of scientists. At some point, however, the questions must arise (and have repeatedly arisen in the literature) of the level of precision with which the Prandtl-vonKármán-Millikan theory (which of course was originally conceived as, and cannot be other than, an approximation) is valid, how quickly it is approached in the limit of Reynolds number tending to infinity, and whether it can be improved with higher-order corrections in some kind of formal asymptotic expansion. It was recently observed [13] that a correction proportional to the pressure gradient, considered as a small perturbation, can explain the difference between circular-pipe, plane-parallel-duct and plane-Couette flow and accord them to a single universal law. Purpose of the present article is to explore how much more can be extracted, from numerical and experimental data available today, in the way of an answer to the other general questions, and how much of this answer can be captured in an analytical interpolation of empirical data.

We must preliminarily remark that an analysis purporting to extract finer and finer corrections purely from empirical data must look at the data at increasing levels of magnification, and will soon meet a point where what is seen is no longer data but just either systematic or random error. One must therefore resist the temptation to overinterpret the data. Doing so necessarily involves subjective judgement, and will force us to make some bold conjectures about what in the available data constitutes the physical phenomenon and what constitutes error. Like all conjectures, some of these may in the future turn out to be wrong, and what is identified here as error may eventually turn out to be an even finer physical feature (or viceversa). I am sure the authors of those data will understand that the term “error” is not intended to convey any negative connotation and is used here in the same sense as in the theory of measurement error.

Nevertheless, as will be clear at the end of this paper, the velocity profiles extracted from today’s available numerical simulations and experiments are affected by a quantifiable error which unfortunately is insufficiently small for some parts of the present analysis. While this detracts nothing from the considerable effort and patience put by the authors in their computations and experiments (an error bar must well be placed somewhere), it does point out the usefulness to pursue further refinements.

In this paper extensive use will be made of velocity profiles taken from the literature, which will be tagged by authors’ initials and Reynolds number. Tags and the references to be credited for each profile are collected in Table 1. We warmly thank all the authors involved for making their precious data files publicly available on the web.

1.1 The classical theory of parallel turbulent flow

The general time-averaged, fully developed velocity profile $u(z)$ in an infinitely long straight duct of a given cross section is a function of the pressure gradient $p_{x}$, density $\rho$, kinematic viscosity $\nu$ and a characteristic height $h$ of the cross-section. Momentum balance links the pressure gradient to the mean wall shear stress $\tau_{w}$ as

[TABLE]

where $D_{H}$, the hydraulic diameter, is 4 times the ratio of area to perimeter of the cross section. It may sometimes be convenient to assume $D_{H}$ itself as the characteristic dimension $h$, but for the time being we shall let them be distinct. $h$ will be chosen to be the distance from the wall to the symmetry axis, i.e. half the distance between walls of a plane duct or the radius of a pipe. Nondimensionalization allows the velocity profile to be expressed as a function of only two variables:

[TABLE]

where $u^{+}=u/u_{\tau}$, $u_{\tau}=\sqrt{\tau_{w}/\rho}$, $z^{+}=zu_{\tau}/\nu$, and $h^{+}=hu_{\tau}/\nu$, quantities denoted by a + being also commonly denoted as measured in “wall units”. The dimensionless height $h^{+}$ in wall units coincides with the shear-based Reynolds number $\mathrm{Re}_{\tau}=hu_{\tau}/\nu$; one or the other name will be more explicative depending on context.

The classical asymptotic theory of the turbulent velocity profile [15, see also 13] distinguishes three regions in the limit of $\mathrm{Re}_{\tau}=h^{+}\rightarrow\infty$:

a “wall” layer, for $z^{+}\ll h^{+}$, where $u^{+}$ is a function of $z^{+}$ only,

[TABLE]

• 2.

a “defect” layer, for $z^{+}\gg 1$, where the velocity defect $u^{+}(h^{+})-u^{+}(z^{+})$ is a function of $Z=z/h=z^{+}/h^{+}$ only,

[TABLE]

• 3.

an “overlap” layer, for $1\ll z^{+}\ll h^{+}$, where both (3) and (4) must simultaneously be valid. According to the standard argument by Millikan [15], their simultaneous validity implies that in the overlap layer the velocity profile is logarithmic,

[TABLE]

where the reciprocal of the $A$ coefficient is known as von Kármán’s constant $\kappa=1/A$.

In addition, equations (3) and (5), operating in a layer much thinner that the characteristic height $h^{+}$, must also be independent of geometry, which is why (3) is named the “universal law of the wall”. By contrast function $F(Z)$ in (4) is allowed to vary from one geometry to another.

The layered structure of the turbulent velocity profile can also be recast as a uniformly valid approximation, by summing (3) and (4) and subtracting their common asymptotic behaviour (5), i.e.

[TABLE]

or equivalently

[TABLE]

$G(Z)$ represents the deviation of velocity in the defect layer from the logarithmic law, defined as

[TABLE]

with $C=\lim_{Z\rightarrow 0}[F(Z)+A\log(Z)]$. By definition of $C$ (not the same $C$ as used in [13]), one has $G(0)=0$. It also follows that $G(1)=C$ [because from (4) $F(1)=0$] and that

[TABLE]

This last relationship is at the foundation of the friction law linking friction coefficient to Reynolds number (e.g., [25] p. 515).

Function $G(Z)$ is known as the “law of the wake” from the classical work of Coles [8], who conjectured on an empirical basis that, in a boundary layer, $G(Z)$ might have a universal shape approximately similar to a mixing layer (or half of a wake), up to a multiplicative constant.

1.2 Pressure-gradient correction to the overlap layer

Several authors have questioned the validity of the logarithmic law, and proposed alternatives, on the basis that (5) provides a somewhat unsatisfactory fit of empirical (both experimental and numerical) data; alternately one may also conclude that (5) is valid but the Reynolds number of present-day direct numerical simulations (or even experiments) is still insufficiently high to see the true asymptotic behaviour and estimate von Kármán’s constant $\kappa$. Such pessimistic deductions are substantiated by plots of the logarithmic derivative of $u^{+}$, a sample of which is given in Figure 1. According to the logarithmic law (5), a plot of the derivative of $u^{+}$ against $\log z^{+}$ should flatten to a constant $A=\kappa^{-1}$ in some intermediate region to be identified with the overlap layer, and admittedly whether it does in Figure 1 can be doubtful.

More insight can be gained if, instead of plotting the velocity profile or its logarithmic derivative for the same geometry at different Reynolds numbers, we plot them at the same Reynolds number for different (circular or plane) geometries. This was done in Figure 1 of our previous article [13], reproduced as Figure 2 here, plotting together three turbulent velocity profiles at the same Reynolds number $\mathrm{Re}_{\tau}=1000$, one for pipe flow from [10], one for plane duct flow from [16] and one for plane Couette flow from [24]. Not only the curves are not straight lines, as is even more clearly brought out in the companion plot of their logarithmic derivative, but in addition the maximum value of velocity is considerably different in each case and so is its average slope. Upon looking at Figures 1-2 one may:

conclude that the logarithmic law is invalid, or not yet valid for $\mathrm{Re}_{\tau}\leq 1000$, except perhaps for Couette flow; 2. 2.

strike a straight line through the velocity plot, assuming that deviations will be somehow compensated, and choose an average slope value as the best approximation of the coefficient $\kappa^{-1}$, to be definite if very crude let us say

[TABLE] 3. 3.

select the minimum at $z^{+}\simeq 70$ common to Figures 1 and 2b, which is also the largest value of $z^{+}$ where all the curves more or less coincide, as the “true” value of $\kappa^{-1}$ and assume that this minimum will keep the same value as the profile becomes flatter and flatter with increasing $\mathrm{Re}_{\tau}$ (a trend made plausible by a comparison of velocity profiles with changing $\mathrm{Re}_{\tau}$ such as in Figure 1).

Route 1 is probably the one that over time elicited the strongest opinions. Considerations about the lack of straightness of logarithmic velocity profiles similar to Figures 1 and 2 have emerged from both experimental data (extending to much larger Reynolds numbers) and numerical simulations. Alternatives to the logarithmic law have been investigated, the most influential and articulate proponent probably being Barenblatt [4].

Route 2 leads to the conclusion (championed by [19, 20], although with a subtler argument and more delicate differences than intentionally exaggerated here) that the velocity profile is indeed logarithmic but not universal. In the relatively low-Reynolds example of Figure 2, (9) would give

$\kappa=0.431$ for plane Couette flow,

• 2.

$\kappa=0.386$ for plane duct flow,

• 3.

$\kappa=0.350$ for circular pipe flow.

Route 3 is tempting, and we gave it some consideration before finding the more satisfactory solution that was described in [13]. It corresponds to $\kappa=0.434$, about the same value as found for Couette flow above.

But in fact, none of these three alternatives turned out to be correct. They are here just to exemplify in how many different ways the same data may sometimes be interpreted. In [13] a two-term asymptotic expansion of the overlap layer was introduced as

[TABLE]

where $A_{0}=\kappa^{-1}$, and $A_{1}$ is a new universal constant quantifying the effect of pressure gradient upon the turbulent velocity profile. The geometry parameter $g$, defined as $-hp_{x}/\tau_{w}=4h/D_{H}$, accounts for the different pressure gradient that corresponds to the same wall shear stress in each geometry, and takes on the values of $g=2$ for circular pipe flow, $g=1$ for pressure-driven flow in a plane duct and $g=0$ for turbulent Couette flow between countermoving plane walls. Equation (10) collapses the three curves of Figure 2 much better than (5), as shown in Figure 2 of [13], and a universal description of the overlap layer is thus restored.

1.3 Joint profile analysis

A digression must be made on the classical method used by experimentalists since Nikuradse [21] to estimate the theoretical parameters of the velocity profile (5). This method is in fact none of the above, but rather to draw together in a single plot the velocity profiles obtained for different Reynolds numbers, and then graphically extract a common trend. The famous Figure 24 of Nikuradse [21], reproduced as Figure 3 here, does so in a graphical format that uses marker symbols to denote experimental points; the unconscious effect is that the overall point cloud strikes the eye much more than the shape of individual velocity profiles. Nikuradse draws two straight lines, one corresponding to $\kappa=0.400$ that more or less cuts through the entire cloud and one corresponding to $\kappa=0.417$ that privileges measurements closer to the wall; the latter he believes to be more accurate because, as he explicitly remarks, his measurements near the centerline must be biased as the profiles do not end with a flat slope there. (We shall return to this observation later.) Nevertheless, the individual profiles are not visually distinguishable; but they become so if the same data are replotted as lines (Figure 4), which is easy to do because Nikuradse wisely provided his data in tabular form.

Had he recognized the importance of this information, he might have discovered the wake function twenty years before Coles.

In fact, from Figure 4 we can observe that each individual curve is consistently slanted upwards with respect to their common trend, so that the slope estimated using a single Reynolds number can be nonnegligibly higher than the one estimated by joining all the data together. Their difference is due to the wake function $G(Z)$, as will be expanded upon later, and should be kept in mind when comparing data of different sources, possibly extracted by different methods. In fact, using more than one Reynolds number at the same time is a fortunate way to compensate for errors induced by $G(Z)$.

On the other hand, figures such as Figure 3 can give the impression that all the point scatter is to be disregarded as measurement error. This way to interpret experimental data may have generated the legend that with the logarithmic law alone “excellent agreement is obtained not only for points near the wall but for the whole range up to the axis of the pipe” [25, p. 509], tantamount to saying, in more recent terminology, that the wake function $G(Z)$ is negligible for pipe flow. Instead, it can be seen from (10) (and will be confirmed by examples in §3) that the wake function of pipe flow is twice as large as the one of plane parallel flow.

1.4 Matched asymptotic expansions

Equation (10) is consistent with a pair of matched asymptotic expansions à la van Dyke [29] of the wall and defect layers, in which an asymptotic expansion in powers of $\mathrm{Re}_{\tau}^{-1}$ of the former is matched to a Taylor series in powers of $Z$ of the latter.

As documented in the reviews [6, 22, 23], several authors, starting from the observation that the ratio between the viscous length scale $l=\nu/u_{\tau}$ and the outer length scale $h$ in Prandtl-vonKármán-Millikan theory is ${h^{+}}^{-1}=\mathrm{Re}_{\tau}^{-1}$, have proposed matched asymptotic expansions in powers of this parameter, modelled after the matched asymptotic expansions of laminar boundary-layer theory111 The choice of $\mathrm{Re}_{\tau}^{-1}$ as the expansion parameter is not obvious, though, because the ratio $u_{\tau}/U$ (which a posteriori turns out to be of the order of $(\log\mathrm{Re}_{\tau})^{-1}$) is also a valid competitor for this role. $u_{\tau}/U$ was assumed as the basic expansion parameter of the boundary layer, for instance, by Mellor [14]. On the other hand, the choice of $\mathrm{Re}_{\tau}^{-1}$ as the expansion parameter is made plausible by the general form of Millikan’s matching condition and by the inclusion of Galilean invariance among its premises.. Typically an inner expansion for the wall layer $0\leq z^{+}\ll h^{+}$ is written as

[TABLE]

and an outer expansion for the defect layer ${h^{+}}^{-1}\ll Z\leq 1$ as

[TABLE]

The derivative of each of the $f_{i}$ functions is then assumed to admit an asymptotic expansion in negative powers of $z^{+}$ for $z^{+}\rightarrow\infty$ (producing a logarithmic term in $f_{0}$ and possibly higher terms) and each of the $F_{i}$ functions an expansion in positive powers of $Z$ for $Z\rightarrow 0$, leading to an order-by-order intermediate matching (e.g. [1, 2, 7, 28, 31]) in the overlap layer.

It can be remarked that the assumption of an algebraic approach of the $f_{i}$’s towards their asymptotic behaviour (a series of negative powers of $z^{+}$ summed to as many positive powers as needed to fit the $F_{i}$’s) is at variance with the classical example of the laminar boundary layer, where this approach is exponentially fast. Whether an algebraic or an exponential behaviour better fits the turbulent flow is an interesting question we shall return to in §5. On the other hand, there is no prominent alternative to a regular expansion in positive powers of $Z$ as a suitable representation of the $F_{i}$’s.

Of course one has to be very careful when deducting an expansion like (11) or (12) from purely empirical data, even more careful than when deriving it from a differential equation, because the finite accuracy of the data will quickly be exceeded with increasing number of terms allowed. For purposes of interpolation, an expression inspired by the uniform approximation (6) may actually be more useful. This can be simply written as

[TABLE]

There is more than one, totally equivalent, possibility of defining the individual terms, because the intermediate overlap behaviour can indifferently be assigned to the $f$ or to the $G$ terms. For instance, the $A_{1}$ term in (10) can be considered as pertaining to $G_{0}$, with the advantage that (6) is unchanged at this order, or to $f_{1}$, doing which emphasizes that the pressure-gradient correction is of first order in the wall layer.

In what follows we shall adopt the convention of associating the logarithmic and constant terms of (10) to $f_{0}$, thus imposing that $G_{0}(0)=0$, and the linear term to $G_{0}$, thus imposing that $f_{1}$ = o$(z^{+})$ for $z^{+}\rightarrow\infty$. Notice that equivalently we might have associated the linear term to $f_{1}$, but in that case (13) to first order would have become

[TABLE]

with $H_{0}(z)$ defined as

[TABLE]

Both $H_{0}(Z)$ and $G_{0}(Z)$ vanish for $Z\rightarrow 0$ but the first does so much faster; either could be apt to be named “law of the wake”, but $G_{0}(Z)$ is more directly comparable to Coles’ definition.

In the following sections we shall endeavour to extract the leading terms of (13) from empirical (numerical and experimental) data, aiming at identifying the wake functions of the three parallel geometries of plane duct, circular pipe and Couette flow, and the universal wall function common to all three. Let us begin by examining each geometry in turn.

2 Plane parallel duct flow

Flow driven by a pressure gradient between two infinite parallel plates (in short, plane-duct flow, or turbulent Poiseuille flow) is the simplest configuration where turbulent DNS (direct numerical simulation) can be performed, and the one for which the majority of DNS data are available, in contrast to flow in a parallel long pipe (turbulent Hagen-Poiseuille flow) for which the majority of experiments are performed. Since, as will be seen, the errors (intended as deviations of one data set from another) of DNS tend to be smaller we shall examine DNS data first, a sample of which is displayed in Figure 5, with the purpose of separating functions $f(z^{+})$ and $G(Z)$ from each other in equation (6). We shall start with the wake function which, being the smaller contribution, can afford a lower relative accuracy and still significantly improve the overall result.

2.1 Estimating the wake function

Examples of the wake function $G(Z)$ for pipes and channels are respectively given in Figs. 3,4 of Panton [23], where they are obtained by difference between each measured, or numerically obtained, actual profile and a preassumed logarithmic law. This is the simplest way to obtain an estimate, and we shall return to it near the end of this paper, but the result is of course sensitive to the a-priori assumption of a log law with specific coefficients as Panton himself notes.

However, if just the uniform approximation (6) is assumed to hold (for the time being, without any higher-order correction), and no log law is explicitly specified, there still is an unbiased way to identify its component functions $f(z^{+})$ and $G(Z)$. These functions can (in the absence of measurement or computation error) be extracted exactly from a pair of experimental or numerical profiles, measured on the same geometry at two different Reynolds numbers. Comparing the results obtained from more than one pair of Reynolds numbers then provides a test on the validity of (6) itself.

Let $u^{+}_{1}$, $u^{+}_{2}$ denote two velocity profiles at Reynolds numbers $\mathrm{Re}_{\tau,1}=h^{+}_{1}$ and $\mathrm{Re}_{\tau,2}=h^{+}_{2}>h^{+}_{1}$. According to (6),

[TABLE]

This is a linear system of two equations in two unknowns from which we can easily eliminate $f(z^{+})$ by subtraction. On defining

[TABLE]

we obtain

[TABLE]

with $m=h^{+}_{1}/h^{+}_{2}$. $\Delta G_{12}(Z)$ itself already constitutes a rough approximation of $G(Z)$ when $h_{2}\gg h_{1}$, since by definition $G(0)=0$ and by continuity $G(mZ)\rightarrow 0$ when $m\rightarrow 0$. But in fact, if (18) is regarded as a functional equation with $G(Z)$ as the unknown, we can just as well solve this functional equation exactly. Recursively applying the substitution $G(Z)\simeq\Delta G_{12}(Z)$ to the second term of (18) provides the tentative solution

[TABLE]

Since, by replacement of the independent variable,

[TABLE]

it is easily verified that (19) is actually the exact solution of (18). The summation converges quickly, and its numerical truncation to some finite $n$ is fairly straightforward. Once $G(Z)$ is calculated, $f(z^{+})$ and its asymptotic, logarithmic or otherwise, behaviour can be immediately obtained (for $0\leq z^{+}\leq h^{+}_{2}$) from (16b) as

[TABLE]

A collection of curves obtained by applying (17,19) to the numerical data of [9, 16, 24] are exhibited in Figure 6. As may be seen, a general common shape of the defect-layer correction $G(Z)$ tends to be recognizable, as expected if (6) is correct, but at the same time a large variability is present, not only between different Reynolds numbers but also between the results of different authors at the same Reynolds number, the data of [24] being the most extreme222Notice that it is the subtraction of two profiles which amplifies the difference. On a standard $u^{+}$-$z^{+}$ plot the results of these same authors look hardly distinguishable. In this sense the present technique might also do double duty as a criterion of numerical convergence.. This is a first indication that the precision (perhaps just the length of time averaging) of the numerical simulations available today still leaves room for improvement; it also hampers the distinction whether the difference between curves in Figure 6 is to be ascribed rather to higher-order corrections of the type (13), or to inaccuracy of the data. Nonetheless one can reckon that, if the data of [24] are excluded, the maximum of $G(Z)$ in Figure 6 is of the order of $0.7$, and arrive at this number without any presumption of a logarithmic law or of the value of von Kármán’s constant. Again with the exception of the data of [24], the vertical range of oscillation of $G(z)$ is of the order of $\pm 0.1$ on top of a velocity $u^{+}$ of the order of $20$, or about $\pm 0.5\%$. This can be easily ascribed to fluctuation of the finite-time numerical data averaging, as it would if (6) were exact, although interestingly such fluctuation remains smooth across $Z$.

Another feature of Figure 6 that may strike one’s attention is that all curves end far from horizontal at $Z=1$, which clashes with Coles’ idea of a half-wake (or mixing-layer) profile. This is the “corner defect” whose discovery is ascribed by Panton [23] to Lewkowicz [12], although according to Musker [18] it was noticed much earlier; with hindsight it is totally to be expected, because the $Z$-derivative of (5) at $Z=1$ is a constant $\kappa^{-1}$ independent of Reynolds number, and therefore $G^{\prime}(1)=-\kappa^{-1}$ is required in order to have the velocity derivative $u_{z}$ vanish on the symmetry axis. The corner defect is of the same order of magnitude as the whole velocity defect, and remains so when $\mathrm{Re}_{\tau}\rightarrow\infty$, thus somewhat invalidating the traditional image of Coles’ law of the wake as a monotonic quarter of a sine wave. (Although, it should be remarked, the latter was conceived with the boundary layer in mind, and Figure 6 might very roughly look like a sine extended over more than a quarter period.)

Despite the relatively small correction provided by $G(Z)$ and the large relative imprecision in its estimate, when $f(z^{+})$ is extracted from one of these curves according to (20), we obtain the third curve in Figure 7, with a very evident progress towards a range of nearly constant logarithmic derivative.

In particular, even if $z^{+}df/dz^{+}$ is not yet as constant as one would desire it to be in order to be able to estimate $\kappa$ precisely, it appears clearly for the first time that the constant range begins at $z^{+}\simeq 200$ (and is distinct from the minimum at $z^{+}\simeq 70$).

The result of Figure 7 should eliminate any doubt that the Reynolds number of present numerical simulations is large enough to spot a logarithmic law; therefore, future research efforts should be directed at improving the accuracy of DNS more than at further increasing the Reynolds number.

Despite the strikingly different appearance of the extrapolated and original logarithmic-derivative profile, an enlargement of Figure 7 shows that a systematic drift is still present. Figure 8 reports the reciprocal of the logarithmic derivative, which can be interpreted as a straight-out estimate of $\kappa$, in a zoomed-in subrange centered about $\kappa=0.4$. As can be seen, the precision of this plot is still insufficient to distinguish what the value of $\kappa$ is, except insofar as it is comprised between, say, $0.37$ and $0.42$. Nevertheless the error is not at random: we see a pattern of seemingly linear drift which in addition appears to change slope, and even sign, with Reynolds number (and, even if not shown in this figure, with data source). A possible interpretation and a way to overcome this drift will be given in next section.

2.2 Eliminating the wake function

Upon closer inspection of Figure 6, one possible pitfall of the previous method comes to the eye: at small $Z$ the curves do not follow outer scaling. While this is not in itself unexpected, being the mark of higher-order corrections of the form (13), it may produce an error in the solution of (19) which propagates from smaller to larger $Z$. In order to circumvent this error we can adopt a complementary approach which moves from large to small $Z$. This will bear some resemblance to the joint profile analysis of §1.3, the one that underlies the data elaboration of many experiments; even if the composite expansion (6) was formulated at a later time than the experimental method became established, it provides in fact a solid background for this method.

If the difference of (16a) and (16b) is taken at equal $Z$, rather than at equal $z^{+}$, instead of eliminating $f(z^{+})$ we obtain the complementary effect of eliminating $G$:

[TABLE]

In the defect layer, if $f\simeq A\log z^{+}+B$ is assumed to hold with an unknown value of $A$, $\Delta f_{21}$ is expected to tend to a constant

[TABLE]

therefore an alternate, easy way to estimate $A$, and von Kármán’s constant $\kappa=A^{-1}$, is to plot the quantity $\Delta f_{21}/\log(h^{+}_{2}/h^{+}_{1})$ as a function of $Z$ and look out for a constant range. A few such plots are displayed in Figure 9 for various pairs of velocity profiles.

Just as was the case in Figure 8, the estimate is consistent with the presence of a range where the plotted quantity is fairly constant (notice that the vertical scale is quite enlarged), but does not exhibit the level of precision required to locate the value of von Kármán’s constant precisely. As shown by the overlaid dashed lines, however, this figure displays a remarkable pattern: once inside the defect layer, all curves tend to align along straight lines that radiate from a common point located at $A_{0}=2.55$, or $\kappa=0.392$.

Despite the straight-line pattern in Figure 9 is very evident, its origin is not so, but it can be remarked that, whereas a single horizontal line would be expected on the basis of (6), (10), and the universality of the $A_{0}$ and $A_{1}$ constants, the observed behaviour becomes predictable if a common value of $A_{0}$ but a different value of $A_{1}$ is used to fit each curve. The best-fit values of $A_{1}$ are listed in Table 2.

A hint at the origin of this variability in $A_{1}$ comes from the observation that there are differences between results of different authors at the same Reynolds number, and the dependence on Reynolds number itself looks somewhat erratic. These are indications that variations in $A_{1}$ are more likely ascribed to numerical error than to the physics of turbulence. Why the error should settle along straight lines is not obvious, but we can offer the following speculation: the slowest decaying mode that perturbs the time averaging is tantamount to a modification in $A_{1}$, or in the interaction between imposed pressure gradient and imposed flow rate. That this or another be the correct interpretation, somewhat astonishingly numerical error produces straight lines in this plot. In what follows we shall proceed based on the conjecture that

Equation (10) is the correct expression for the overlap layer, our present best approximation of $\kappa$ is $0.392$, and different numerical simulations err in their estimate of $A_{1}$ according to Table 2.

This conjecture is tantamount to an, admittedly bold, estimate of the numerical error in order to subtract it, and the reader may or may not want to accept its implications, but some intriguing deductions can be drawn as will soon appear. Well aware that only a set of more accurate simulations will be able to clarify the matter in the future, we shall present the results of such deductions both with a single common value of $A_{1}$ and with the values of $A_{1}$ taken from Table 2.

Whatever the origin of the pattern in Figure 9, the value of $A_{0}=1/0.392$ is bound to play a special role, and as will be seen all our further tests in this paper are consistent with it being our present best approximation of $\kappa$. It can also be noted that the estimates of $A_{1}$ in Table 2 tend to cluster about unity, and most importantly, the pair of profiles that in Figure 9 turn up to have a difference with almost flat slope, i.e. the best adherence to a logarithmic law (HJ550 and HJ2E3), have in Table 2 a value of $A_{1}$ very close to unity. All this leads us to further infer that the correct value may be exactly $A_{1}=1$, as will also be confirmed by the comparison with cylindrical pipe flow in §3. The values of $\kappa=0.392$, $A_{1}=1$ and $B=4.48$, giving rise to the logarithmic law

[TABLE]

were reported in [13] as our present best estimates.

An immediate consequence of allowing a different value of $A_{1}$ for each DNS is that a very precise identification of the corner-defect function $H(Z)$ of (15) ensues. This is brought out in Figure 10, where (10) is subtracted from each profile with $A_{0}=1/0.392$ and $A_{1}$ taken from Table 2. The resulting estimates of $H(Z)$ are nearly identical (compare the accord in Figure 10 to the imprecision in Figure 6), and for the purpose of interpolation can be fit by the power-law

[TABLE]

It follows from (15) that, within this approximation, Coles’ wake function $G(Z)$ can be usefully expressed as

[TABLE]

Once again, this function has a different shape than traditionally has been associated with Coles’ law of the wake (a quarter sine wave with flat slope both at $Z=0$ and at $Z=1$, or a similarly shaped cubic polynomial). It must be remarked, however, that Coles [8] drew this shape with the boundary layer in mind, which might behave differently from plane-duct flow in this respect.

Eventually, when the complete wake function (24) is subtracted from the velocity profile so as to extract the wall function $f(z^{+})$, the outcome is Figure 11 with $A_{1}=1$ or Figure 12 with $A_{1}$ taken from Table 2, in either case a very substantial improvement over Figure 5 as far as agreement to a logarithmic law is concerned.

3 Pipe flow

Numerical velocity profiles for circular pipe flow are less abundant and extend to a smaller Reynolds number than for flow in a plane duct; in exchange the majority of experimental data pertain to this geometry.

3.1 DNS data

Figure 13 shows the unmanipulated velocity profiles of El Khouri et al. [10] and Wu and Moin [30] up to $Re_{\tau}\gtrsim 1000$. As can be seen, they depart even more pronouncedly from a straight line than the plane-duct profiles did in Figure 1.

When the wake-extraction technique of §2.1 is applied to these profiles, the result is Figure 14.

Here the variability is larger than in Figure 6, mostly due to the lower range of Reynolds number but also probably to a shorter period of time averaging; nonetheless a very similar shape of the wake function can be spotted, except that the vertical scale is doubled. As was explained in [13], this is no accident and can be correlated to the fact that the pressure gradient is also doubled in this geometry. Here just as in Figure 6, an initial kink arises on the wall scale; this kink requires a word of caution since it might be responsible, through error propagation in (19), of a possibly artificial vertical shift of the whole curve.

When the wake-elimination technique of §2.2 is applied to pipe flow DNS results, the result is Figure 15. The pattern of straight lines looks forced upon this figure, again an indication of the lower range of Reynolds numbers but also probably of a generally larger residual numerical error, and it cannot honestly be spotted without a comparison with Figure 9. Nonetheless, nontrivially this plot is compatible with a straight-line pattern emanating from the same value of $\kappa$ as the plane profiles. This becomes fairly evident if, instead of taking differences between profile pairs, we accept as an ansatz to be proved that the log law coincides with the one found for plane flow, including its coefficients $\kappa=0.392$ and $B=4.48$, and, in Figure 16, subtract from each profile the log law (22) and the corner-defect function.

There should now be hardly any doubt that the difference approaches a straight line with a slope close to $gA_{1}=2$ ($g$ being the geometry factor defined in §1.2, $g=2$ for pipe flow). The best-fit slopes are reported in Table 3.

As to the corner-defect function,

this is shown in Figure 17 and can be approximated as

[TABLE]

When eventually the complete wake function

[TABLE]

is subtracted from each velocity profile, the estimate of the universal law of the wall in Figure 18 ensues, which should be compared to the untreated Figure 13 as far as adherence to a logarithmic law is concerned.

3.2 Experimental data

Experiments on turbulent pipe flow extend up to a Reynolds number $\mathrm{Re}_{\tau}\simeq 10^{5}$, which is a hundred times larger than the largest value yet achieved in DNS. On the other hand they present generally larger variability, and the higher the Reynolds number the worse they resolve the wall region for obvious instrument-size constraints. Lack of the wall region renders the method of §2.1 inapplicable, since (19) needs a reasonably good interpolation there.

The wake-elimination method of §2.2 does not suffer from this limitation, since profile differences are only taken in the defect layer, but error amplification by such differences renders the result, visible in Figure 19, of little quantitative value (even though the barycenter of the whole figure is not far from the by now expected $A_{0}=2.55$).

On the other hand, if again the same wake function (26) is subtracted as obtained from previous cases, Nikuradse’s original data in Figure 4 become those in Figure 20,

and exhibit a close fit to the the logarithmic law (22) with coefficients $\kappa=0.392$ and $B=4.48$; a fit which is even more convincing insofar as these coefficients were extrapolated from a Reynolds number a hundred times smaller. In fact, at closer inspection the individual velocity profiles in Figure 20 are now no longer slanted and much more parallel to the log law, the main discrepancy being a deviation downwards in the rightmost part of each curve corresponding to the corner defect. Nikuradse himself comments that his measurements are not completely reliable near the center.

Actually Nikuradse [21] concluded that his most reliable estimate was $\kappa=0.417$, but he narrowly missed a more accurate result with his Figure 18, where he looks for the dependence of velocity on Reynolds number at a constant position rather than the opposite. Had he plotted $u/u_{\tau}$ as a function of $\mathrm{Re}_{\tau}$ for constant $Z$, instead of $u/U$ as a function of $\mathrm{Re}_{\text{bulk}}$ for constant $Z$, he would have implicitly eliminated the contribution of the wake function $G(Z)$ and obtained the present Figure 21, where the best fit slope corresponds to $\kappa=0.3904$.

Similar considerations apply to the Superpipe experiment of [26], whose raw velocity profiles (actually, “raw” after the fine corrections described in [26]) are reported in Figure 22a. The upward slant is even more evident than in Nikuradse’s.

When the wake function (26) is subtracted, these become the profiles displayed in Figure 22b, with again a sizeable progress towards the same logarithmic law (22). The main discrepancy is this time a slight residual upwards deviation in the rightmost part of each curve, which when compared to Nikuradse’s excessive downwards deviation corroborates the idea that such residual deviations, in either case smaller than deviations of the original data, might actually be measurement errors.

Taken all together, the different experimental curves tend to satisfactorily tighten up onto the logarithmic law (22) when the wake function, including the pressure-gradient correction, is subtraced. It must barely be reminded that this is not a logarithmic law fitted to these particular data, but one having the same numerical coefficients (both $A_{0}$ and $B$) as in all our other plots.

4 Couette flow

Turbulent Couette flow has been much less intensely investigated than pipe or plane-duct flow, and we have the DNS of [24] as our only source. Nevertheless it is an essential contribution to our comparison, as it completes the picture with a third value of the geometry factor $g$. ($g=2$ for pipe flow, $g=1$ for plane-duct flow, $g=0$ for turbulent Couette flow.)

The wake-estimation curves for Couette flow, shown in Figure 23,

exhibit a behaviour markedly different from previous cases, with a much lower positive range ($+0.1$ instead of $+0.7$, effectively zero in the face of foreseeable data errors). The lack of a significant positive overshoot, in the light of the explanation given in [13], can be immediately associated with the value of $g$ being zero. In the corner-defect region we see negative values, consistent with a negative corner-defect function $H(Z)$ which in this case is not offset by any pressure-gradient effect. It is also quite evident that $\mathrm{Re}_{\tau}=260$ is likely to be too low, and only the difference between the $\mathrm{Re}_{\tau}=500$ and $\mathrm{Re}_{\tau}=1000$ profiles can be trusted here.

The difference between these two profiles in the wake-elimination plot of §2.2 is shown in Figure 24.

This is very satisfactorily consistent with $A_{0}=1/0.392$ and $gA_{1}=0$, even if such values could not have been evinced with precision from this plot alone.

The comparison between Couette velocity profiles (untouched, since their pressure-gradient correction is zero) and the by now standard logarithmic law (22) is shown in Figure 25.

Finally, the corner-defect (actually coincident with the wake) function is shown in Figure 26, as extracted by subtracting the log law from each Couette-flow profile.

As can be seen, the initial kink that appeared in Figure 23 really was an artifact of the method; according to Figure 26 the wake function is negative and very nearly monotonic, becomes significant already for $Z\geq 0.4$ (earlier than in the other two geometries), and for practical purposes it is very well approximated by the smoothed-ramp function

[TABLE]

Whereas the general shape of the corner-defect function was similar for plane and pipe flow, for Couette flow we see a significant difference. This is unsurprising since here the profile must adapt to a different (odd rather than even) symmetry. It is nonetheless a remarkable coincidence that the maximum negative value of the corner-defect correction ($\sim-0.5$), occurring at the centerline $Z=1$, is similar in all three geometries.

5 The law of the wall

Let us now more closely examine the wall layer, with the purpose of obtaining a description, and possibly some analytical interpolating formula, of the universal law of the wall $f_{0}(z^{+})$ and of its first-order correction $f_{1}(z^{+})$ in (13). Particularly interesting will be how quickly the sum of these functions approaches its overlap behaviour (10), so as to know in which range of $z^{+}$ and with what accuracy the simpler equation (10) can be relied upon. In doing so we shall zoom in on finer and finer details of the velocity profile, and it should be kept in mind that error of whatever origin is correspondingly amplified and we are venturing in a path where physical effect and numerical (or measurement) error may become indistinguishable, something which is inevitably bound to happen sooner or later when empirical data are being analysed at deeper and deeper resolution.

Once the logarithmic law (22) is taken for granted, the difference between each velocity profile and this log law can be highlighted by rescaling plots similar to Figures 10,17,26 in wall units rather than in outer units. Such a plot for plane parallel flow is shown in Figure 27a, where only the logarithmic law is subtracted without any pressure-gradient correction,

and in Figure 27b where also the pressure-gradient correction is included as

[TABLE]

Several noticeable features appear. First of all the agreement between profiles obtained by different authors, at different Reynolds numbers, and in different geometries is very good on the scale of Figure 27a and even better in Figure 27b, thus showing that there really is a universal law of the wall. In Figure 27a we see a linear (as it would appear had we not adopted a logarithmic scale) divergence after $z^{+}\simeq 100$; as can by now be imagined, this linear divergence is caused by the pressure gradient, and vanishes when the pressure-gradient correction is additionally subtracted in Figure 27b. The residual uncertainty in Figure 27b is of the order of $\pm 0.02$ in its left half and $\pm 0.1$ its right half, a region which will later require some caution.

Secondly, the wall-function difference exhibits a well defined minimum at $z^{+}\simeq 2.6$ and a well defined maximum at $z^{+}\simeq 42$, both of which very convincingly scale with wall units. The overshoot of the wall function before approaching the logarithmic law is a remarkable feature, and was highlighted as such before [17], since it eludes the instinctive expectation of a smooth, one-sided approach. It also explains the generalized slight overestimation of $\kappa$ (underestimation of $A_{0}$) in the past, as people were prone to apply a logarithmic law already for $z^{+}\gtrsim 30$. Nikuradse’s original preference for his estimate of $0.417$ over $0.4$ in Figure 3 may perhaps also be ascribed to this tendency.

One may be tempted to interpret this overshoot as the presence of a “mesolayer”, conjectured in [3] to possess an intermediate length scale $\sqrt{h\nu/u_{\tau}}$; however such an interpretation, which may perhaps be left open by a figure such as 27a, looks unwarranted in the light of Figure 27b where the position and shape of the maximum clearly stay fixed in wall units with changing Reynolds number.

5.1 Extrapolation to $\mathrm{Re}_{\tau}=\infty$

Whereas it should be apparent that any one of the profiles included in Figure 27b (or an analytical interpolation of one of those) can already yield a good approximation of the leading-order law of the wall $f_{0}(z^{+})$, or of its increment $\Delta f_{0}(z^{+})$ to be added to (10), trying to extract the first-order correction $f_{1}$ drives us to the next level of magnification.

Truncating (11) to its first two terms

[TABLE]

shows that, for constant $z^{+}$, velocity $u^{+}$ is expected to be a linear function of $\mathrm{Re}_{\tau}^{-1}$ (up to some value of $z^{+}$ where (29) ceases to be valid). An empirical verification whether this is so is provided by plots such as Figures 28a–c, which display velocity as a function of $\mathrm{Re}_{\tau}^{-1}$ for three different values of $z^{+}$. As can be seen, the individual points reasonably fit a straight line but not perfectly, and in fact the difference between results of different authors becomes visible on this scale even at positions relatively near to the wall.

Though Figure 28 does not afford a precise evaluation of the first-order correction to the wall function, it does allow us to extract its order of magnitude and to state that its absolute value, of $0.1$ in the worst case for $\mathrm{Re}_{\tau}=500$, is small enough not to be practically very relevant. On the other hand, the linear extrapolation to $1/\mathrm{Re}_{\tau}=0$ of the data in Figure 28 looks like the appropriate way, in the absence of additional random errors, to extract the most accurate approximation of $f_{0}(z^{+})$ from the available data files. This extrapolation is shown, in a neighbourhood of its maximum, in Figure 29 together with a straightforward arithmetic mean of the curves in Figure 27b. It should be remarked that extrapolation may be dangerous, as it amplifies the non-deterministic part of the error. The arithmetic mean, on the contrary, damps random data errors but introduces a systematic bias.

As can be seen, the difference between the two estimates crosses zero near their maximum, located at $z^{+}\simeq 42$, and grows larger afterwards where, however, the spread among data in Figure 27b was also larger. Caution suggests as a compromise to adopt the extrapolated curve up to the crossing point and the averaged curve thereafter, and this will be our present best approximation of the law-of-the-wall increment $\Delta f_{0}$ defined by (28).

5.2 Interpolation

In lack of a closed physical model of mean turbulent flow, obvious practical convenience has always made analytical interpolations of the law of the wall desirable. A very early such interpolation was devised by Spalding [27], in the form of a monotonic function satisfying the end conditions $f_{0}(0)=0$, $f^{\prime}_{0}(0)=1$, $f^{\prime\prime}_{0}(0)=0$, $f^{\prime\prime\prime}_{0}(0)=0$, and $f(z^{+})\rightarrow\kappa^{-1}\log(z^{+})+B$ for $z^{+}\rightarrow\infty$. He obtained this result by first inverting the logarithmic law (5) as $z^{+}=\operatorname{e}^{\kappa(u^{+}-B)}$, and then subtracting as many terms of the Taylor series of this exponential as required to zero the appropriate number of derivatives:

[TABLE]

Spalding’s formula is elegant but implicit, as (30) needs to be inverted to give $u^{+}$ as a function of $z^{+}$. As will be seen below, it is also the least accurate.

An explicit formula was provided by Musker [18], who started with a rational (ratio of two polynomials) expression of Reynolds stress, designed to satisfy appropriate asymptotic conditions at $z^{+}=0$ and $z^{+}\rightarrow\infty$, and then analytically integrated it to extract the velocity profile as the formula

[TABLE]

Neither Spalding nor Musker used other data than behaviour near zero and infinity to constrain their interpolation. More recently, Monkewitz et al. [17] adopted a rational Padé representation of the Reynolds stress formally similar to Musker’s, but using a larger number of free coefficients, and fitted it to a set of experimental profiles measured in a zero-pressure-gradient boundary layer. (Like Musker’s, this was part of a more extensive effort aimed at describing both the wall and wake regions of a turbulent boundary layer. Nonetheless, the part of their analysis concerning the wall layer also applies to parallel flow insofar as the wall layer is universal.) Monkewitz et al. [17] split the Reynolds stress in the sum of two contributions, using a Padé function P23, ratio of a polynomial of degree $2$ and one of degree $3$, for one and a Padé function P25 for the other, and determined the coefficients of these two rational functions from a combination of boundary conditions and experimental data points. An explicit expression of the velocity profile requires analytically integrating the Reynolds stress, and takes the form of a sum of several logarithms and arctangents of polynomials (Eq. (6) of [17]).

In order to develop an interpolation of the data in Figure 29, we preliminarily observe that the integral of a rational function will always have an asymptotic behaviour at infinity given by the integral of its Laurent power series, namely one that can be expressed as a single logarithm of a first-degree polynomial plus another rational function. Therefore, in order to set up a generic formula with free coefficients apt to fit empirical data, we would get a more easily manageable expression if, instead of integrating a rational function like in the above examples, we wrote $f_{0}(z^{+})$ as the sum of a logarithm and a Padè approximant and fitted this sum to the data.

As an additional consideration, either way of using rational functions still restricts their asymptotic behaviour to an algebraic form, which is another name for a Laurent series. This is a very specific kind of asymptotic behaviour, involving analyticity in the complex plane, and among the differential equations of mathematical physics it is a prerogative of the Laplace equation in open space. In many other physical phenomena, such as those governed by constant-coefficient differential equations (including the Laplace equation between walls) an exponential approach to their behaviour at infinity is more common. Particularly in the matched asymptotic expansion of boundary layers [29], to which the asymptotic extensions of Millikan’s theory are inspired, approach of the inner-layer solution to its own asymptotic behavior is exponential.

Whereas a Padé interpolation over an infinite interval will eventually converge given enough free coefficients, it will do so slowly if it does not match the appropriate asymptotic behavior. In order to empirically decide whether the data in Figure 29

look more like exponential rather than algebraic, a simple device is to replot them on a logarithmic vertical scale, as done in Figure 30. The appearance of about a decade of nearly constant slope provides an indication that the asymptotic decay of $\Delta f_{0}$ is more likely to be exponential than algebraic. (Or at least, aside from any theoretical speculations, an exponential interpolant is more likely to succeed in interpolating this curve.) Balancing this consideration with the undoubtful practical convenience of rational functions, our final choice is a P22 Padé approximant multiplied by an exponential, summed to the logarithm of a first-degree polynomial:

[TABLE]

Imposing the initial conditions $f_{0}(0)=0$ and $f^{\prime}_{0}(0)=1$, together with the asymptotic values $\kappa=0.392$ and $B=4.48$, and performing a least-square fit of the remaining 5 coefficients to the data of Figure 29, yields the values listed in Table 4.

A plot of (32) is shown in Figure 31, compared to the composite extrapolated and averaged data of Figure 29 and to the past approximations of [17, 18, 27], each represented as a wall-function increment $\Delta f_{0}$ obtained by subtracting one and the same log law (22). The error of the interpolant (32) with coefficients from Table 4, undistinguishable from zero in Figure 31, is further displayed in Figure 32 and is of the order of $\pm 0.01$ in wall units.

6 Roundup and conclusion

After going through the discussion in the previous sections, one sees that our present best approximation to the mean velocity profile in different geometries of parallel flow still fits within Coles’ shrewd formula

[TABLE]

where, however, there are now explicit interpolations available of the universal wall function $f_{0}(z^{+})$ and of each specific wake function $G(Z)$. $f_{0}(z^{+})$ is for all practical purposes well represented by its analytical interpolant (32), and $G(Z)$ is given by

[TABLE]

$gZ$ being the effect of the pressure gradient333More precisely, as shown in [13], the effect of the pressure gradient must be written as $A_{1}gZ$ where $A_{1}$ is a universal constant, derived from the same kind of dimensional argument as von Kármán’s constant is, the best approximation of which is presently $A_{1}=1$., with the geometry factor $g$ equal to $2$ for pipe flow, $1$ for plane duct flow, and [math] for turbulent Couette flow, and $H(Z)$ being the corner-defect function approximated by, respectively, (25), (23), and (27) for each geometry.

An even better representation must in principle be the higher-order uniformly-valid formula

[TABLE]

but an estimate of $f_{1}(z^{+})$ sits on the edge of what can be gotten within the accuracy of present simulations and experiments. This additional contribution to $u^{+}$ is, at any rate, of less than $0.1$ absolute value for $\mathrm{Re}_{\tau}\geq 500$.

In order to extract the wall function $f_{0}(z^{+})$ from empirical data we have, in practice, first extracted an approximation of the wake function $G(Z)$ from the comparison of more than one velocity profile at different Reynolds numbers, and then subtracted $G(Z)$ from one of those. The technique to do this was explained in §2.1, consisting in the solution of a functional equation involving the difference of two velocity profiles taken at the same wall-units coordinate $z^{+}$. The method is computationally easy, and provides a visual image of the wake function that does not require any assumption other than the validity of (33); it also shows that the extracted wall function actually becomes logarithmic for $z^{+}\gtrsim 200$, thus the range of Reynolds numbers for which numerical simulations are available is at all sufficient in order to try and estimate the value of von Kármán’s constant.

Unfortunately, this estimate still contains too much error (data uncertainty) to provide a sufficiently precise value of $\kappa$. Therefore in §2.2 we proposed a different method that, by taking differences of velocity profiles at the same outer coordinate $Z$ rather than at constant $z^{+}$, mathematically eliminates the wake function even if doing so hides its shape. The outcome of this method, Figure 9 in the case of plane-duct flow, is expected to be a set of superposed horizontal straight lines pointing at the true value of $\kappa$, but much to our surprise it turned out to be close to a pencil of inclined straight lines radiating from a single point located at $A_{0}=1/0.392$. The inclination of such lines suggests that the different profiles have somewhat different wake functions, and in the light of the theory [13] which associates the slope of the wake function with the pressure gradient, it is as though the numerical error mimicked somewhat different pressure gradients. (Alternately, this might be a true physical effect that makes (33) invalid, but the difference being of the same order of magnitude as the difference between results of different authors makes us lean towards numerical error as the source.)

Pending a verification with more accurate data in the future, we settled on $\kappa=0.392$ and $A_{1}=1$ as our present best estimates. We then proceeded to examining various datasets of plane-duct, circular-pipe and turbulent-Couette flow in the light of this logarithmic law and extracted both the wake function of each geometry and a common law of the wall, for which the compact analytical interpolation (32) was provided.

Having decided on what $f_{0}(z^{+})$ should look like, we can now go full turn and test the validity of our previous assumptions, and at the same time the accuracy of the data, by reconstructing the wake function of each velocity profile as $G(Z)=u^{+}-f_{0}(z^{+})$. This is done in Figures 33-36, where the theoretical wake function, as deduced in the corresponding section of this article, is also superposed. As can be seen there is a considerable residual spread of the order of $\pm 0.2$ for a plane duct and even larger for pipe flow (we cannot comment about spread for Couette flow since a single data set is available in this case), but the fact that the predicted wake function, using one and the same wall function for all three geometries, fits within the data cloud in all cases is a nontrivial confirmation.

The open question is, of course, is the residual spread just caused by error (and therefore will decrease with future improvements in numerical and experimental measurement techniques), or does it accomodate some yet undiscovered physical feature? Closer inspection reveals a possible pattern in these figures, behind the evident random scatter, of general decrease of $G$ with Reynolds number. This pattern might point at a slightly revised estimate of $\kappa$, or at a different shape of the wake function (perhaps a quadratic at small $Z$, as Coles [8] seemed to expect, rather than a linear function?), or at the presence of higher-order corrections of the form (13). But it might also point at some systematic error in the numerical techniques (non-uniform discretization or size of the periodic box for instance, just to mention two whose influence changes with Reynolds number) or in the experimental techniques. It is in the human nature to look for patterns, and strong is the temptation to speculate further, but such a chase for finer and finer details can never end and easily go astray. We prefer to content ourselves with the present level of investigation and wait for future data improvements, with the only final remark that the techniques and plots described in this paper can also become a useful tool to analyse the data themselves and their accuracy during their retrieval.