A recently developed mixing length model of the turbulent shearing stress in wall bounded flows has been used to formulate a universal velocity profile (UVP) that provides an effective replacement for the widely used Coles wall-wake formulation. Comparisons with both direct numerical simulation and experimental data demonstrate the ability of the profile to approximate a wide variety of wall-bounded flows. The UVP is uniformly valid from the wall to the boundary layer edge and for all Reynolds numbers from zero to infinity. There is no presumption of logarithmic dependence of the velocity profile outside the viscous wall layer so the profile can accurately approximate low Reynolds number turbulent boundary layers. The effect of a pressure gradient is included in the UVP through the introduction of a modified Clauser parameter that correlates well with the parameters that determine the wake portion of the velocity profile. The inherent dependence of the UVP on Reynolds number, extended to include the effect of pressure gradient, enables it to be used as the basis of a new method for integrating the von Kármán boundary layer integral equation for a wide variety of attached wall bounded flows. To illustrate its application, the UVP is used to determine the zero-lift drag coefficient of the Joukowsky 0012 and NACA (National Advisory Committee for Aeronautics) 0012 airfoils over a wide range of chord Reynolds numbers.

In volume III of the Durand series on Aerodynamic Theory, Prandtl1 points out that the boundary layer approximation derived on a flat wall can be extended to a curved surface. He writes: “In this case, a coordinate system can be chosen in which the arc length along the fixed surface can be introduced as the x-coordinate, and the perpendicular distance from the surface as the y-coordinate.” The idea is illustrated in Fig. 1.

FIG. 1.

Curved surface boundary layer coordinate notation. The wing shown is a Joukowsky 0012.

FIG. 1.

Curved surface boundary layer coordinate notation. The wing shown is a Joukowsky 0012.

Close modal

Using this construct, a number of simplified approaches were developed to produce approximate solutions suitable for estimating the drag of planar aerodynamic shapes. For laminar boundary layers, the most well-known integral methods are due to Thwaites2 and Pohlhausen.3 The latter was extended recently by Majdalani and Xuan.4 There are also a limited number of methods that integrate the Kármán equation for turbulent boundary layers. One of the most widely used is the method due to Head.5 He assumes the existence of two universal functions, F and G, for the turbulent boundary layer. The first function relates a modified boundary layer shape factor to the conventional shape factor, and the second relates the conventional shape factor to the entrainment velocity at the outer edge of the boundary layer. When these functions are combined with the friction law of Ludwieg and Tillman,6 all of the basic characteristics of the boundary layer can be determined. The main weakness is that the functions F and G are determined from experimental correlations that exhibit very limited universality and the Ludwieg–Tillman friction law is valid only over a narrow range of Reynolds numbers.

Low Mach number turbulent flow over a flat plate at high Reynolds number is described by the incompressible boundary layer approximation (1)

(1)

subject to the no-slip and free stream conditions

(2)

where u, ue, and v are mean velocities.

In the region near the stagnation point, the boundary layer approximation breaks down (Fig. 2). Outside of this region, if the surface radius of curvature is large compared to the boundary layer thickness, no significant pressure variation will occur in the wall-normal direction and wall curvature plays no role in the boundary layer approximation.7 Theoretical underpinning of Prandtl's idea can be found in the invariance of the boundary layer approximation, (1), under the Lie translation group with an arbitrary function of x added to the y coordinate (Cantwell,8 Sec. 10.5, p. 301).

FIG. 2.

Coordinate notation near the wing leading edge.

FIG. 2.

Coordinate notation near the wing leading edge.

Close modal

The free stream velocity $ue(x)$ is determined by solving Laplace's equation for the potential flow about the body with the Kutta condition applied at the trailing edge. With $ue(x)$ known, along with a suitable model of the turbulent shear stress $u′v′¯$, the calculation of boundary layer characteristics can proceed.

Boundary layer integral methods all use the von Kármán9 equation (3) derived by integrating the boundary layer equations (1) over the height of the boundary layer

(3)

The friction velocity is defined as

(4)

The displacement thickness is defined as

(5)

and the momentum thickness is

(6)

where δh indicated in Fig. 1 is a boundary layer thickness defined in Cantwell10 as the equivalent channel half height thickness. The thickness, δh, will be discussed further in Sec. II B.

Wing drag is calculated using the construction in Fig. 3. The differential force acting in the $x̃$ direction on, say, the wing upper surface is

(7)
FIG. 3.

Drag force on a differential section of the wing surface.

FIG. 3.

Drag force on a differential section of the wing surface.

Close modal

The viscous drag coefficient based on the chord of a symmetric wing of the type shown in Fig. 1 is

(8)

where the local (based on ue) friction coefficient is

(9)

The universal velocity profile comes from the solution of the fully developed channel (or pipe) flow equation

(10)

where the pressure gradient is a constant, Cantwell.11

Express Eq. (10) in wall normalized coordinates

(11)

Integrating (10) with respect to y+ from the lower wall of the channel, the result is the stress balance

(12)

where $Rτ$ is the friction Reynolds number

(13)

The turbulent shear stress is modeled using classical mixing length theory (von Kármán;12 see also Prandtl13 and van Driest14). Let

(14)

When (14) is substituted into (12), the result is a quadratic equation for $du+/dy+$. We use the positive root

(15)

The singularity at the wall where the mixing length, $λ$, goes to zero can be easily removed as follows.15,16 Multiply (15) by one,

(16)

Carry out the multiplication in (16) and canceling terms, the result, equivalent to (15), but with the singularity removed is

(17)

where $du+dy+(0)=1$ and $du+dy+(Rτ)=0$. The limiting velocity gradients at the wall and centerline are

(18)

independent of the choice of $λ(y+)$.

Equation (17) is integrated from the wall to y+ to obtain the velocity profile in the form of an integral dependent on the non-dimensional mixing length function $λ(y+)$ at a given $Rτ$

(19)

where $u+(0)=0$ and $u+(Rτ)=ue/uτ$. At low Reynolds number, $λ→0$ and Eq. (19) approaches the laminar pipe/channel flow solution

(20)

where, in the laminar limit, $Rτ=(2ueδh/ν)1/2$.

The mixing length model introduced by Cantwell11 to approximate pipe data is

(21)

The model (21) contains five free parameters. The constant k is closely related to the Kármán constant. The parameter a constitutes a wall damping length scale. The wall model is similar to the exponential decay proposed by van Driest14 except for the exponent m that helps to determine the shape and thickness of the near wall velocity profile. Near the wall

(22)

The outer flow term in the denominator of (21) includes a length scale b proportional to the fraction of the boundary layer thickness where wake-like behavior begins as well as an exponent n that helps shape the outer part of the profile.

The reason for making a distinction between k and the Kármán constant is that the UVP accurately approximates low Reynolds number turbulent velocity profiles that have no distinct logarithmic layer and a fairly complicated dependence on k. Whereas at high Reynolds number ($Rτ>2000/k$),11k becomes simply a scale factor on u+, y+ and the damping length scale, a. In that regime, a distinct logarithmic dependence in the velocity profile begins to emerge and k can be viewed as equivalent to the usual definition of the Kármán constant, κ.

The UVP provides an effective replacement for the classical wall-wake formulation.17 The velocity profile (19) with the mixing length model (21) is uniformly valid over $0≤y≤δh$ and $0≤Rτ<∞$. Therefore, there is no need for a buffer layer function and there is no discontinuity in the velocity derivative at the outer edge of the boundary layer. At very low Reynolds number, the velocity profile reverts to the laminar solution. The profile is directly connected to a model of the turbulent shear stress that can be used in computations based on the full Reynolds Averaged Navier–Stokes equations. Finally, the UVP and the shear stress model can be used to study turbulence characteristics such as turbulent kinetic energy (TKE) production.

The classical wall-wake formulation is silent as to the Reynolds number at which the formulation actually applies. Whereas for $Rτ≥2000/k$, the integral form of the UVP, Eq. (19) approaches the explicit form10,11

(23)

Importantly, the shape function, $ϕ$, is independent of $Rτ$. As long as $Rτ$ is above about $2000/k$, the explicit form of the UVP can be used to replace the integral form (19) outside the viscous sublayer and buffer layer, $y+>132$ (see Table II in Subrahmanyam et al.17).

The shape function for the zero pressure gradient (ZPG) boundary layer using the parameter values in Table I is shown in Fig. 4. A best fit eighth-order polynomial is overlaid on $ϕ$. In Sec. IV, the shape function will be extended to include pressure gradient effects and the eighth-order fit will be extended to a bivariate polynomial. In practice, the shape function is determined by simply evaluating the UVP at any $Rτ≥2000/k$, say $Rτ=105$ or 106, to produce $ϕ$. Any velocity profile with $Rτ≥2000/k$ can then be generated by integrating Eq. (19) only out to about $y+=132$.17 The rest of the profile out to $y+=Rτ$ is generated using Eq. (23). Complete details are given in Cantwell10,11 and Subrahmanyam et al.17

TABLE I.

Average model parameters with standard deviation for basic wall flows. Ranges of $Rτ$ for each flow are as follows: Pipe (3327–530 023), Channel (550–8016), ZPG boundary layer (1343–17 207).

Flow$k¯$σk$a¯$σa$m¯$σm$b¯$σb$n¯$σn
Pipe (21 profiles) 0.4092 0.0057 20.0950 0.381 1.6210 0.0379 0.3195 0.0157 1.6190 0.1204
Channel (7 profiles) 0.4086 0.0179 22.8673 1.599 1.2569 0.0292 0.4649 0.0485 1.3972 0.1213
ZPG boundary layer (11 profiles) 0.4233 0.0068 24.9583 0.663 1.1473 0.0373 0.1752 0.0060 2.1707 0.2238
Flow$k¯$σk$a¯$σa$m¯$σm$b¯$σb$n¯$σn
Pipe (21 profiles) 0.4092 0.0057 20.0950 0.381 1.6210 0.0379 0.3195 0.0157 1.6190 0.1204
Channel (7 profiles) 0.4086 0.0179 22.8673 1.599 1.2569 0.0292 0.4649 0.0485 1.3972 0.1213
ZPG boundary layer (11 profiles) 0.4233 0.0068 24.9583 0.663 1.1473 0.0373 0.1752 0.0060 2.1707 0.2238
FIG. 4.

The universal velocity profile shape function using the parameter values in Table I evaluated at $Rτ=106$ overlaid by an eighth-order polynomial fit (dashed line).

FIG. 4.

The universal velocity profile shape function using the parameter values in Table I evaluated at $Rτ=106$ overlaid by an eighth-order polynomial fit (dashed line).

Close modal

In this way, the UVP can be evaluated up to arbitrary large Reynolds numbers enabling the structure of wall bounded flows to be explored in the limit of infinite Reynolds number.10,11,15,18

The high Reynolds number friction law is generated by evaluating (23) at the boundary layer edge

(24)

shown as the dashed line in Fig. 5. The blue curve in Fig. 5 is the exact friction law: Eq. (19) evaluated at $y+=Rτ$. Note that the two curves merge around $Rτ=2000/k$.

FIG. 5.

The universal velocity profile friction law for a zero pressure gradient boundary layer. The dashed line is the high Reynolds number limit generated by evaluating Eq. (23) at $y+=Rτ, y/δh=1$.

FIG. 5.

The universal velocity profile friction law for a zero pressure gradient boundary layer. The dashed line is the high Reynolds number limit generated by evaluating Eq. (23) at $y+=Rτ, y/δh=1$.

Close modal

The integral method introduced here relies on the UVP to define the velocity profile over the surface in question. In Subrahmanyam et al.,17,19 the UVP is shown to provide an accurate approximation to zero pressure gradient (ZPG) boundary layer computational data from Simens et al.,20 Borrell et al.,21 Sillero et al.,22 and Eitel-Amor et al.23 The UVP accurately approximates the ZPG experimental data of Klebanoff,24 DeGraaff and Eaton,25 and Baidya et al.26 as well as the adverse pressure gradient data of Perry and Marusic.27 In the present paper, the UVP will be shown to accurately approximate the favorable pressure gradient (sink flow) experimental data of Jones.28

In fact, as noted earlier, the velocity profile that accurately approximates all these flows is fundamentally the same channel flow solution. One flow is distinguished from another only by the values of the five empirical constants $(k,a,m,b,n)$ that characterize it. Since $∂u/∂y=0$ at the edge of the boundary layer, the UVP defines an unambiguous overall thickness for the boundary layer; the equivalent channel half height, δh, referred to earlier. For a given flow geometry, the Reynolds number dependence of the model parameters $(k,a,m,b,n)$ tends to be quite weak and so average values of the model parameters can provide a good approximation to the velocity profile over a wide range of Reynolds numbers. In Cantwell,10 a boundary layer integral method utilizing the universal velocity profile was used to generate the integral properties of the zero pressure gradient boundary layer on a wall over the full Reynolds number range $0≤Rτ<∞$. In the present paper, this method is generalized to boundary layer flows with pressure gradients and the parameters of the mixing length model used to approximate the wake portion of the velocity profile in the UVP are allowed to vary depending on the pressure gradient. Importantly, Subrahmanyam et al.17 showed that the wall parameters (k, a, m) exhibit very little change in the presence of an adverse pressure gradient.

To facilitate the later discussion, four functions will be defined here; F0, F1, F2, and F3. As noted above, the boundary layer friction law is generated by evaluating the UVP at $y+=Rτ$

(25)

This is the blue curve in Fig. 5 just discussed. The local displacement thickness Reynolds number expressed in wall units is

(26)

and the local momentum thickness Reynolds number is written as

(27)

For the UVP integral method, we will also need the derivative of the momentum thickness Reynolds number

(28)

The model parameters $(k,a,m,b,n)$ for the mixing length function (21) are selected by minimizing the sum of total squared error between a given data profile and the universal velocity profile (19) using the cost function

(29)

Parameter values $(k¯,a¯,m¯,b¯,n¯)$ averaged over a wide range of Reynolds numbers for smooth wall channel, pipe and zero-pressure-gradient boundary layer flow are shown in Table I along with the variances of each parameter over the related dataset. The corresponding ranges of $Rτ$ are given in the caption. Throughout the present paper, the ZPG boundary layer model constants are assumed to be the boundary layer average values shown in Table I, $(k¯,a¯,m¯,b¯,n¯)=(0.4233,24.9583,1.1473,0.1752,2.1707)$ from Table I in Subrahmanyam et al.17 Generally speaking the parameters $(k,a,m,b,n)$ depend, at most weakly on the Reynolds number for a given flow geometry and at high Reynolds number the dependence appears to be especially weak. The evidence for this is the analysis of $(k,a,m,b,n)$ for the Princeton Superpipe data11 that spans three orders of magnitude in the Reynolds number.29 Table I includes the variance in the optimal parameters over the set of profiles for each flow giving a sense of the variation over the ranges of $Rτ$ for each flow geometry. For the zero pressure gradient boundary layer, the variation in k, a, m, and b is less than $4%$ while the variation of n is about $10%$.

The statistical variation in parameter values indicated in Table II implies a similar variation in the wall friction derived from the UVP. Figures 6 and 7 provide a relatively simple measure of the extent of the variation in friction associated with variations in the parameter values by simply showing the friction law with parameter values increased by one standard deviation and decreased by one standard deviation with no attempt to account for any correlation that might exist between parameters. The variation in friction coefficient shown in Fig. 7 translates directly to a comparable variation in calculated drag. For later reference, at $Rτ=5000$, the friction coefficient based on the mean parameters in Table I is $Cf=0.002 378$. With the mean parameters all increased by one standard deviation $Cf+=0.002 463$ and with the parameters decreased by one standard deviation $Cf−=0.002 277$ for an overall variation of approximately $7.5%$. In Sec. VIII where airfoil drag comparisons are discussed, this uncertainty in Cf will be discussed again.

TABLE II.

Run data, Reynolds number, optimal model parameters, and RMS error for adverse pressure gradient boundary layer datasets from Perry and Marusic.27 Initial free stream values are $u∞=10$ and $u∞=30$ m/s. Channel half height thickness for these data are at $u=0.998ue$; $δh=δ0.998$. Reprinted with permission from Subrahmanyam et al., J. Fluid Mech. 933, A16 (2022). Copyright 2022 Author(s), licensed under a Creative Commons Attribution (CC BY) License.17

x (m)$ue(ms)$$Rδ1$$Rδ2$ββc$δ998(m)$$Rδ998$$Rτ$$(ueuτ)$kambn$u+rms$
1.20 10.361 3 165 2 282 0.0 0.0 0.031 79 21 439 912 23.51 0.4287 25.18 1.1528 0.212 2 2.111 0.120
1.80 9.976 5 226 3 734 0.65 1.115 0.050 19 32 606 1285 25.37 0.4239 24.59 1.1583 0.183 9 1.705 0.253
2.24 9.256 6 410 4 342 1.45 2.432 0.055 43 33 456 1195 28.00 0.4223 24.70 1.1460 0.123 7 2.461 0.163
2.64 8.588 8 606 5 517 2.90 4.760 0.070 55 39 406 1252 31.47 0.4255 24.79 1.1121 0.093 31 2.399 0.152
2.88 8.155 11 235 6 879 4.48 7.223 0.086 34 46 043 1337 34.44 0.4296 25.54 1.1173 0.073 83 3.178 0.220
3.08 7.896 12 397 7 213 7.16 11.326 0.092 63 47 598 1248 38.13 0.4302 25.20 1.0938 0.061 51 4.578 0.231
1.20 30.704 8 772 6 564 0.0 0.0 0.033 53 64 807 2461 26.34 0.4095 23.21 1.1161 0.228 5 1.743 0.0705
1.80 29.054 12 401 9 073 0.71 1.230 0.044 15 80 849 2870 28.17 0.4088 23.25 1.1468 0.172 2 1.922 0.0895
2.24 27.035 16 307 11 587 1.39 2.378 0.055 26 94 275 3137 30.05 0.4112 23.42 1.1710 0.133 1 2.335 0.0942
2.64 25.150 21 634 14 736 2.74 4.606 0.069 68 110 700 3373 32.82 0.4035 22.87 1.1352 0.106 6 2.663 0.0984
2.88 23.885 25 854 17 020 3.96 6.567 0.080 54 121 760 3471 35.08 0.4018 22.80 1.1213 0.090 7 3.164 0.1183
3.08 22.908 31 767 20 052 6.07 9.901 0.093 73 136 290 3587 37.99 0.4083 23.53 1.1339 0.074 2 4.069 0.1673
x (m)$ue(ms)$$Rδ1$$Rδ2$ββc$δ998(m)$$Rδ998$$Rτ$$(ueuτ)$kambn$u+rms$
1.20 10.361 3 165 2 282 0.0 0.0 0.031 79 21 439 912 23.51 0.4287 25.18 1.1528 0.212 2 2.111 0.120
1.80 9.976 5 226 3 734 0.65 1.115 0.050 19 32 606 1285 25.37 0.4239 24.59 1.1583 0.183 9 1.705 0.253
2.24 9.256 6 410 4 342 1.45 2.432 0.055 43 33 456 1195 28.00 0.4223 24.70 1.1460 0.123 7 2.461 0.163
2.64 8.588 8 606 5 517 2.90 4.760 0.070 55 39 406 1252 31.47 0.4255 24.79 1.1121 0.093 31 2.399 0.152
2.88 8.155 11 235 6 879 4.48 7.223 0.086 34 46 043 1337 34.44 0.4296 25.54 1.1173 0.073 83 3.178 0.220
3.08 7.896 12 397 7 213 7.16 11.326 0.092 63 47 598 1248 38.13 0.4302 25.20 1.0938 0.061 51 4.578 0.231
1.20 30.704 8 772 6 564 0.0 0.0 0.033 53 64 807 2461 26.34 0.4095 23.21 1.1161 0.228 5 1.743 0.0705
1.80 29.054 12 401 9 073 0.71 1.230 0.044 15 80 849 2870 28.17 0.4088 23.25 1.1468 0.172 2 1.922 0.0895
2.24 27.035 16 307 11 587 1.39 2.378 0.055 26 94 275 3137 30.05 0.4112 23.42 1.1710 0.133 1 2.335 0.0942
2.64 25.150 21 634 14 736 2.74 4.606 0.069 68 110 700 3373 32.82 0.4035 22.87 1.1352 0.106 6 2.663 0.0984
2.88 23.885 25 854 17 020 3.96 6.567 0.080 54 121 760 3471 35.08 0.4018 22.80 1.1213 0.090 7 3.164 0.1183
3.08 22.908 31 767 20 052 6.07 9.901 0.093 73 136 290 3587 37.99 0.4083 23.53 1.1339 0.074 2 4.069 0.1673
FIG. 6.

The blue line is the UVP friction law for a zero pressure gradient boundary layer using mean parameters $(k¯,a¯,m¯,b¯,n¯)$ from Table II. The dashed line is the UVP friction law with each parameter increased by one standard deviation, $(k¯+σk,a¯+σa,m¯+σm,b¯+σb,n¯+σn)$. The dotted line is the UVP friction law with each parameter decreased by one standard deviation, $(k¯−σk,a¯−σa,m¯−σm,b¯−σb,n¯−σn)$.

FIG. 6.

The blue line is the UVP friction law for a zero pressure gradient boundary layer using mean parameters $(k¯,a¯,m¯,b¯,n¯)$ from Table II. The dashed line is the UVP friction law with each parameter increased by one standard deviation, $(k¯+σk,a¯+σa,m¯+σm,b¯+σb,n¯+σn)$. The dotted line is the UVP friction law with each parameter decreased by one standard deviation, $(k¯−σk,a¯−σa,m¯−σm,b¯−σb,n¯−σn)$.

Close modal
FIG. 7.

Same as Fig. 6 but expressed in terms of the friction coefficient, $Cf=2/(ue/uτ)2$.

FIG. 7.

Same as Fig. 6 but expressed in terms of the friction coefficient, $Cf=2/(ue/uτ)2$.

Close modal

It should be noted that Fig. 5 and Table II can be compared to Fig. 6 and Table I in Cantwell.10 The average values of $(k,a,m,b,n)$ in Cantwell10 differ slightly from the ZPG parameter values used in the present paper but are the same as in the earlier conference paper by Subrahmanyam et al.19 The reason is that boundary layer data from Baidya et al.26 is included in the average values quoted in Subrahmanyam et al.17 but was not used in determining average optimal parameters in the 2021 papers where it was decided to use direct numerical simulation (DNS) data only, since it reaches all the way to the wall. In retrospect that condition is not necessary, although experimental data do tend to have higher uncertainty. As new data become available, the average optimal parameters will continue to be updated although changes can be expected to be relatively small as time goes on and the sample size grows.

In the presence of a pressure gradient, the parameters $(b,n)$ will be allowed to vary while (k, a, m) will be held fixed at the average values in Table I. The key question of how $(b,n)$ are related to the pressure gradient will be addressed in Sec. IV.

#### 1. The equivalent channel half height thickness

In Fig. 1, the boundary layer thickness, δh, is identified as the equivalent channel half height thickness. The idea behind δh originates in the fact that the UVP accurately approximates not only pipe flow and channel flow but also boundary layer flows with a zero, positive and negative pressure gradient.10,11,17

The UVP has a well-defined outer edge, where $u/ue=1$ and $∂u/∂y=0$. However, these properties of the boundary layer profile are only approached asymptotically. The usual practice is to choose an arbitrary point where $u/ue$ is slightly less than one and $∂u/∂y$ is slightly greater than zero. Probably the most common choice is the boundary layer thickness corresponding to $ue/Ufreestream=0.99$, where Ufreestream is the free stream velocity reported with the data for a given profile.

The UVP can be used to remove some of this arbitrariness. When optimal parameters for the UVP are determined for boundary layer data, the choice of cutoff point for the boundary layer edge has an impact on the rms error. If δh is too large, data points approaching the free stream where the velocity is virtually constant will be included in the profile. The optimization procedure will try to fit these points and this will tend to reduce the accuracy over the whole profile leading to increased error. If δh is too small, the data will be cutoff too far short of the boundary layer edge. The condition $∂u/∂y=0$ will be imposed where the derivative is not quite zero. This also leads to increased error. There is a choice of δh that minimizes the error between the UVP and the data. See pp. 12 and 13 and Figs. 17–23 in Cantwell10 where δh is first defined.

At present, techniques for measuring and reporting data vary widely among researchers and so the boundary layer edge velocity corresponding to δh also varies. Generally, the best fit occurs at $ue=(0.993−0.998)Ufreestream$.

Figure 8 shows how the shape of the universal velocity profile becomes fuller as $Rτ$ is increased by 7 orders of magnitude. At the lowest Reynolds numbers, the profile is essentially the laminar channel/pipe velocity profile. At the highest Reynolds number, large increases produce very small changes in the velocity profile. At extreme Reynolds numbers, the velocity profile approaches plug flow but astronomically large values are needed to approach this state.11

FIG. 8.

The universal velocity profile at $Rτ=1,10,100,1000,104,105,106$, and $107$. Parameters $(k,a,m,b,n)$ are the ZPG boundary layer values from Table I.

FIG. 8.

The universal velocity profile at $Rτ=1,10,100,1000,104,105,106$, and $107$. Parameters $(k,a,m,b,n)$ are the ZPG boundary layer values from Table I.

Close modal

Figure 5 shows the UVP friction law along with low Reynolds number and high Reynolds number limits for a zero pressure gradient boundary layer. Based on Fig. 5, the profile appears to begin becoming fuller at about $Rτ=30$ corresponding to $Rδ1=122, Rδ2=51$ and $Rx=$15 700 with a shape factor H = 2.39. The profile approaches a fully turbulent shape at about $Rτ=500$ corresponding to $Rδ1=2030,Rδ2=1373$ and $Rx=645 000$ with a shape factor H = 1.48. These changes in the profile inferred from Fig. 5 approximately correspond to the first and third profiles in Fig. 8. The evolution of the profile encompasses the nominal transition range from the critical Reynolds number of the Blasius profile, $Rτ=121, Rδ1=520$, and $Rx=104 000$,7 to a nominal transition Reynolds number, $Rδ1=1600, Rτ=392$ and $Rx=474 000$. In other words, the UVP begins to evolve at a Reynolds number well below the instability limit.

The critical Reynolds number for a smooth plate is highly dependent on the free stream turbulence level usually denoted Tu. The question of transition in the UVP boils down to determining what transition assumption one is making when the UVP is used throughout the flow from a laminar leading edge to a fully turbulent boundary layer. To explore this question, the UVP friction law is plotted again in Fig. 9 along with transition data from Schubauer and Klebanoff30 and Coupland.31 Note that large $ue/uτ$ implies low Cf. Extremely low values of Tu can lead to critical Reynolds numbers exceeding $Rτ=2000$ corresponding to approximately $Rex=2 700 000$ while the UVP implies turbulent values of the friction coefficient beginning at about $Rτ≈200$ corresponding to $Rex≈200 000$.

FIG. 9.

The universal velocity profile friction law for a zero pressure gradient boundary layer along with transitional friction data from Schubauer and Klebanoff30 (natural transition $Tu=0.03%$—black) and Coupland31 (ERCOFTAC case T3A $Tu=0.9%$—red, case T3A $Tu=3.0%$—magenta).

FIG. 9.

The universal velocity profile friction law for a zero pressure gradient boundary layer along with transitional friction data from Schubauer and Klebanoff30 (natural transition $Tu=0.03%$—black) and Coupland31 (ERCOFTAC case T3A $Tu=0.9%$—red, case T3A $Tu=3.0%$—magenta).

Close modal

Methods for modeling the transition region on a smooth airfoil, for example, Kaynak et al.,32 available in the aerodynamic simulation suite SU2,33 correlate the transition distance with Tu. In order to accurately determine the drag of a given smooth airfoil, the free stream turbulence needs to be well characterized. In the present paper, in the absence of such data, we will assume the UVP holds beginning at the forward stagnation point while recognizing that that assumes that the critical Reynolds number is close to the minimum required to sustain turbulent flow, $Rτtransition≈200$. In other words, the UVP corresponds to a tripped boundary layer.

Figure 9 can be compared to Fig. 14 in Cantwell11 where the UVP friction law for pipe flow is compared with data over 7 orders of magnitude in the pipe Reynolds number. Based on the discussion here, one can conclude that that comparison is for a UVP pipe profile corresponding to a disturbed pipe entry flow.

In view of the accurate approximation of the UVP to DNS and experimental data, it should be expected that the UVP generates an accurate friction law. This is crucial if the UVP integral method presented here is to be useful and it needs to be confirmed.

Figure 10 shows the UVP friction law from Figs. 5 and 9. Above $Rτ=2000/k$, the UVP friction law is

(30)
FIG. 10.

Five friction laws compared to the UVP (blue line). The orange line is the Ludwieg–Tillman law6 used in Head's method. The magenta line is the friction law developed by Nash.34 The red and black lines are two versions of the Coles–Fernholz law with $(k,C)=(0.41,5.0)$ (red line) and $(k,C)=(0.384,4.1)$ (black line). The brown line is the Spalart–Allmaras law, Eq. (34), deduced from data in Polewski and Cizmas.36

FIG. 10.

Five friction laws compared to the UVP (blue line). The orange line is the Ludwieg–Tillman law6 used in Head's method. The magenta line is the friction law developed by Nash.34 The red and black lines are two versions of the Coles–Fernholz law with $(k,C)=(0.41,5.0)$ (red line) and $(k,C)=(0.384,4.1)$ (black line). The brown line is the Spalart–Allmaras law, Eq. (34), deduced from data in Polewski and Cizmas.36

Close modal

In Fig. 10, the UVP law is compared with five other laws. The Ludwieg–Tillman relation6 used in Head's integral method is

(31)

where H is the shape factor, $H=δ1/δ2$.

Also shown in Fig. 10 is the relation developed by Nash.34 Nash based his friction law on a study of several classical approaches including Ludwieg–Tillman. Nash's approach involves the following three equations:

(32)

Nash's law agrees quite closely with the UVP over a wide range of Reynolds numbers.

Also included in Fig. 10 are two versions of what Nagib et al.35 call the Coles–Fernholz law.

(33)

Finally, Polewski and Cizmas36 carried out a detailed comparison of several turbulence models. In their Fig. 9(a), the skin friction coefficient on a flat plate generated by three codes using the Spalart–Allmaras37 model is plotted up to $Rx=107$. The curves agree very closely. The Spalart–Allmaras friction law at this Reynolds number deduced from this plot is

(34)

Note that for the comparisons in Figs. 10 and 11, the $Rδ1$ and $Rδ2$ that appear in Eqs. (31)–(33) are related to $Rτ$ using the integrals (5) and (6) expressed in terms of $Rτ$, Eqs. (26) and (27). See also the explicit high Reynolds number relations, Eqs. (40) and (43), in Cantwell.10

FIG. 11.

Four friction laws compared to the UVP (blue line) and plotted as Cf in linear coordinates. The magenta line is the friction law developed by Nash.34 The red and black lines are the two versions of the Coles–Fernholz law. The brown line is the Spalart–Allmaras model deduced from data in Polewski and Cizmas.36

FIG. 11.

Four friction laws compared to the UVP (blue line) and plotted as Cf in linear coordinates. The magenta line is the friction law developed by Nash.34 The red and black lines are the two versions of the Coles–Fernholz law. The brown line is the Spalart–Allmaras model deduced from data in Polewski and Cizmas.36

Close modal

The closest to the UVP law in Fig. 10 are the Coles–Fernholz law with Coles' recommended values $(k,C)=(0.41,5.0)$ and the Spalart–Allmaras relation. Outside of the range $1000, the Ludwieg–Tillman law leads to unrealistically high $ue/uτ$ (low Cf) values and will not be considered further.

Figure 11 shows the UVP, Nash, Spalart–Allmaras, and the two Coles–Fernholz laws in linear coordinates over the range $0 which covers the vast bulk of available boundary layer data ($0 corresponding to $0 corresponding to $0). All five curves in Fig. 10 agree within about 5% over the range above $Rτ>1000$. The closest agreement with the UVP law in this Reynolds number range are the laws due to Nash34 and Spalart and Allmaras.37 Note that only the UVP merges with the laminar limit, $Rτ→0$ as shown in Fig. 5.

Nagib et al.35 carried out measurements of the skin friction on a flat plate in the Reynolds number range $10 000. They used that data, along with older data corrected to form a common experimental dataset, in a comprehensive study of several widely used classical friction laws. The result is shown in Fig. 12 where their Fig. 7 is overlaid by the friction law generated by the UVP (blue circles) and the friction law of Nash34 (magenta dashed line). Note that the $Rτ$ and $Rδ2$ domains in Figs. 11 and 12 represent almost identical Reynolds number ranges.

FIG. 12.

Friction coefficient from the UVP compared with Fig. 7 from Nagib et al.35 The friction law of Nash34 is also included (magenta dashed line).

FIG. 12.

Friction coefficient from the UVP compared with Fig. 7 from Nagib et al.35 The friction law of Nash34 is also included (magenta dashed line).

Close modal

Figures 10–12 show that the UVP friction law is consistent with the vast body of friction data in the literature.

Define the overall flow Reynolds number in terms of the free stream velocity $u∞$ and an appropriate body length scale which we will take to be the airfoil leading edge radius of curvature r shown in Fig. 2

(35)

The main drag results are expressed in terms of the chord Reynolds number, Rchord, obtained by multiplying Eq. (35) by $c/r=91.3737$ for the J0012 airfoil or $c/r=63.5870$ for the NACA0012 airfoil. Introduce the length scale r and let $ξ=x/r$. The dimensionless potential flow velocity at the surface of the body is

(36)

Now re-express the Kármán equation (3) in terms of $Rτ$. First, express the momentum thickness, in terms of the momentum thickness Reynolds number (27),

(37)

Differentiate (37) with respect to x

(38)

Using (38) in (3), the Kármán equation is now written

(39)

Normalize x by the nose radius, r, and use (26) and (27) to express δ1 and δ2 in (39) in terms of the displacement and momentum thickness Reynolds numbers. Use the definition (25) for the friction law and the definition (36) for the inviscid velocity.

(40)

The form of the Kármán equation that we need to solve numerically is

(41)

where (28) has been used. The zero pressure gradient form of (41) was used by Cantwell10 to determine integral measures of the ZPG boundary layer over $0≤Rτ<∞$.

Subrahmanyam et al.17 used the UVP to approximate the adverse pressure gradient data of Perry and Marusic27 for the two upstream velocity cases, $u∞=10$ and $u∞=30$ m/s. The near wall parameters $(k,a,m)$ did not vary significantly with the Reynolds number, whereas the outer flow wake parameters (b, n) varied quite a lot. The results from Subrahmanyam et al.17 are shown in Table II along with the fit of the UVP to the Ref. 27 data in Figs. 13 and 14.

FIG. 13.

Reprinted with permission from Subrahmanyam et al., J. Fluid Mech. 933, A16 (2022). Copyright 2022 Author(s) licensed under a Creative Commons Attribution (CC BY) License.17 Comparison between the universal velocity profile and the $u∞=10$ m/s adverse pressure gradient data of Perry and Marusic27 (open red circles).

FIG. 13.

Reprinted with permission from Subrahmanyam et al., J. Fluid Mech. 933, A16 (2022). Copyright 2022 Author(s) licensed under a Creative Commons Attribution (CC BY) License.17 Comparison between the universal velocity profile and the $u∞=10$ m/s adverse pressure gradient data of Perry and Marusic27 (open red circles).

Close modal
FIG. 14.

Reprinted with permission from Subrahmanyam et al. J. Fluid Mech. 933, A16 (2022). Copyright 2022 Author(s) licensed under a Creative Commons Attribution (CC BY) License.17 Comparison between the universal velocity profile and the $u∞=30$ m/s adverse pressure gradient data of Perry and Marusic27 (open red circles).

FIG. 14.

Reprinted with permission from Subrahmanyam et al. J. Fluid Mech. 933, A16 (2022). Copyright 2022 Author(s) licensed under a Creative Commons Attribution (CC BY) License.17 Comparison between the universal velocity profile and the $u∞=30$ m/s adverse pressure gradient data of Perry and Marusic27 (open red circles).

Close modal

The sink flow data from Jones et al.28 will be used to determine how the wake parameters of the UVP, (b, n) vary for favorable pressure gradients. These results are provided in Table III and the comparison between the UVP and the data is shown in Figs. 15–17. The quantity K in Table III is essentially the inverse of the sink flow Reynolds number,

(42)

where $u∞$ is the converging channel entry velocity, x is the coordinate along the lower wall and L is the effective length of the sink flow channel. The optimal value of the UVP parameter k for these data tends to be relatively high compared to ZPG and adverse pressure gradient data, particularly near the entrance to the converging channel and at the lowest Reynolds number. This may be a consequence of the relatively low values of $Rτ$ that characterize the data. Comparable values of k were observed at low $Rτ$ in the pipe data (Cantwell11 Fig. 9, PSP cases 1–5). The mean parameter values in Table I for pipe flow were averaged over PSP cases 6–26.

TABLE III.

Run data, Reynolds number, optimal model parameters, and RMS error for favorable pressure gradient boundary layer datasets from Jones et al.28 Converging channel entry velocities are $u0=5.0,7.5$ and 10.0 m/s. Channel half height thickness for these data are at $u=0.996ue$; $δh=δ0.996$. Kinematic viscosity is $ν=1.51×10−5m/s2$. Note, according to Jones et al.,28$ue/u0=1/(1−x/L)$ and the calibrated sink length is $L=5.60m$.

$K×107$x (mm)$u0(ms)$$Rδ1$$Rδ2$$−β$$−βc$$Rδ996$$Rτ$$(ueuτ)$kambn$u+rms$
5.39 800 5.0 1112 780 0.2436 0.4145 8656 429 20.16 0.4686 26.88 1.2136 0.3931 1.3295 0.098
5.39 1600 5.0 1629 1192 0.3780 0.6548 14 130 681 20.75 0.4630 25.58 1.2408 0.5284 1.1644 0.068
5.39 2200 5.0 1648 1209 0.3880 0.6726 14 806 709 20.90 0.4670 27.89 1.0869 0.5806 1.1343 0.105
5.39 2800 5.0 1946 1449 0.4713 0.8222 18 244 861 21.20 0.4555 26.89 1.1196 0.6934 1.0908 0.089
5.39 3280 5.0 2138 1606 0.5276 0.9241 20 549 960 21.40 0.4445 26.04 1.1563 0.6831 1.1818 0.067
5.39 3580 5.0 2226 1687 0.5500 0.9671 22 468 1049 21.41 0.4448 25.96 1.1920 0.8189 1.1514 0.068
3.59 800 7.5 1496 1069 0.3613 0.4125 11 752 555 21.17 0.4602 26.62 1.1906 0.3738 1.2986 0.092
3.59 1600 7.5 2000 1470 0.3382 0.5868 17 325 798 21.70 0.4644 28.91 1.0413 0.4647 1.1745 0.085
3.59 2200 7.5 2350 1755 0.4030 0.7040 21 366 977 21.86 0.4557 27.30 1.0822 0.5752 1.0988 0.060
3.59 2800 7.5 2827 2155 0.4850 0.8748 27 899 1262 22.12 0.4431 26.19 1.1527 0.7463 1.0917 0.065
3.59 3280 7.5 2928 2245 0.5202 0.9191 29 762 1338 22.25 0.4374 25.93 1.1726 0.7271 1.1931 0.055
3.59 3580 7.5 3027 2339 0.5344 0.9473 32 204 1452 22.17 0.4403 25.78 1.2037 0.8700 1.1579 0.062
2.70 800 10.0 1862 1343 0.3161 0.4092 15 001 690 21.74 0.4564 28.35 1.0045 0.4324 1.1316 0.097
2.70 1600 10.0 2555 1904 0.3421 0.5969 22 546 1013 22.27 0.4529 27.56 1.0538 0.5090 1.1120 0.074
2.70 2200 10.0 2873 2160 0.3908 0.6846 26 114 1164 22.44 0.4525 26.90 1.0669 0.5491 1.0864 0.047
2.70 2800 10.0 3372 2577 0.4689 0.8273 32 769 1444 22.70 0.4402 25.79 1.1307 0.6833 1.0766 0.047
2.70 3280 10.0 3725 2851 0.5297 0.9351 35 899 1564 22.95 0.4364 25.45 1.1227 0.6653 1.0704 0.073
2.70 3580 10.0 3936 3044 0.5575 0.9888 39 990 1746 22.91 0.4368 25.31 1.1607 0.7504 1.0923 0.095
$K×107$x (mm)$u0(ms)$$Rδ1$$Rδ2$$−β$$−βc$$Rδ996$$Rτ$$(ueuτ)$kambn$u+rms$
5.39 800 5.0 1112 780 0.2436 0.4145 8656 429 20.16 0.4686 26.88 1.2136 0.3931 1.3295 0.098
5.39 1600 5.0 1629 1192 0.3780 0.6548 14 130 681 20.75 0.4630 25.58 1.2408 0.5284 1.1644 0.068
5.39 2200 5.0 1648 1209 0.3880 0.6726 14 806 709 20.90 0.4670 27.89 1.0869 0.5806 1.1343 0.105
5.39 2800 5.0 1946 1449 0.4713 0.8222 18 244 861 21.20 0.4555 26.89 1.1196 0.6934 1.0908 0.089
5.39 3280 5.0 2138 1606 0.5276 0.9241 20 549 960 21.40 0.4445 26.04 1.1563 0.6831 1.1818 0.067
5.39 3580 5.0 2226 1687 0.5500 0.9671 22 468 1049 21.41 0.4448 25.96 1.1920 0.8189 1.1514 0.068
3.59 800 7.5 1496 1069 0.3613 0.4125 11 752 555 21.17 0.4602 26.62 1.1906 0.3738 1.2986 0.092
3.59 1600 7.5 2000 1470 0.3382 0.5868 17 325 798 21.70 0.4644 28.91 1.0413 0.4647 1.1745 0.085
3.59 2200 7.5 2350 1755 0.4030 0.7040 21 366 977 21.86 0.4557 27.30 1.0822 0.5752 1.0988 0.060
3.59 2800 7.5 2827 2155 0.4850 0.8748 27 899 1262 22.12 0.4431 26.19 1.1527 0.7463 1.0917 0.065
3.59 3280 7.5 2928 2245 0.5202 0.9191 29 762 1338 22.25 0.4374 25.93 1.1726 0.7271 1.1931 0.055
3.59 3580 7.5 3027 2339 0.5344 0.9473 32 204 1452 22.17 0.4403 25.78 1.2037 0.8700 1.1579 0.062
2.70 800 10.0 1862 1343 0.3161 0.4092 15 001 690 21.74 0.4564 28.35 1.0045 0.4324 1.1316 0.097
2.70 1600 10.0 2555 1904 0.3421 0.5969 22 546 1013 22.27 0.4529 27.56 1.0538 0.5090 1.1120 0.074
2.70 2200 10.0 2873 2160 0.3908 0.6846 26 114 1164 22.44 0.4525 26.90 1.0669 0.5491 1.0864 0.047
2.70 2800 10.0 3372 2577 0.4689 0.8273 32 769 1444 22.70 0.4402 25.79 1.1307 0.6833 1.0766 0.047
2.70 3280 10.0 3725 2851 0.5297 0.9351 35 899 1564 22.95 0.4364 25.45 1.1227 0.6653 1.0704 0.073
2.70 3580 10.0 3936 3044 0.5575 0.9888 39 990 1746 22.91 0.4368 25.31 1.1607 0.7504 1.0923 0.095
FIG. 15.

Comparison between the universal velocity profile and the $K=5.39×10−7$ sink flow data of Jones et al.28 (open red circles).

FIG. 15.

Comparison between the universal velocity profile and the $K=5.39×10−7$ sink flow data of Jones et al.28 (open red circles).

Close modal
FIG. 16.

Comparison between the universal velocity profile and the $K=3.59×10−7$ sink flow data of Jones et al.28 (open red circles).

FIG. 16.

Comparison between the universal velocity profile and the $K=3.59×10−7$ sink flow data of Jones et al.28 (open red circles).

Close modal
FIG. 17.

Comparison between the universal velocity profile and the $K=2.70×10−7$ sink flow data of Jones et al.28 (open red circles).

FIG. 17.

Comparison between the universal velocity profile and the $K=2.70×10−7$ sink flow data of Jones et al.28 (open red circles).

Close modal

The pressure gradient term in (41) can be arrived at from another direction. The strength of the pressure gradient in a turbulent boundary layer is often expressed in terms of the Clauser38 parameter $β=(δ1/τw)(dp/dx)$. See, for example, Perry and Marusic.27 Here, we will introduce a modified Clauser parameter

(43)

The Bernoulli equation is used to express βc in terms of the local free stream velocity, and the thicknesses are expressed in terms of the local thickness Reynolds numbers. The result is identical to the pressure gradient term in the Kármán equation (41)

(44)

The UVP parameters n and b are plotted in Figs. 18 and 19 against βc. The solid line

(45)

approximates the data for $n(βc)$ in Fig. 18 reasonably well with root mean square (rms) error 0.197. Two approximations to the data for $b(βc)$ are shown in Fig. 19. The dashed line is the simpler of the two and is given by

(46)

with rms error 0.0440. Equation (46) gives a good approximation to most of the data but is not accurate near $βc=0$. This is a problem because over most of the airfoil surface the pressure gradient is quite small as can be seen in the lower plot in Fig. 22. This means that βc is also quite small over a substantial part of the airfoil surface. The solid line is given by

(47)

with rms error 0.052 5. The exponential term provides the correction needed to better approximate the data in the neighborhood of $βc=0$ without significantly increasing the error away from $βc=0$.

FIG. 18.

Correlation between the model parameter n and βc for the data in Tables II and III.

FIG. 18.

Correlation between the model parameter n and βc for the data in Tables II and III.

Close modal
FIG. 19.

Correlation between the model parameter b and βc for the data in Tables II and III. The solid line provides a more accurate approximation near $βc=0$.

FIG. 19.

Correlation between the model parameter b and βc for the data in Tables II and III. The solid line provides a more accurate approximation near $βc=0$.

Close modal

The effect of this change can be seen by evaluating Cf at $βc=0$ for $Rτ=104$. If the parameters b and n from Table I, $(b,n)=(0.1752,2.1707)$, are used $Cf=0.002 13$. If Eqs. (45) and (46) are used, $(b,n)=(0.3050,1.4194)$ and $Cf=0.002 38$, which is $12%$ higher than the Table I value and will give almost a $12%$ higher airfoil drag. If Eqs. (45) and (47) are used, $(b,n)=(0.2223,1.4194)$ and the friction coefficient is $Cf=0.002 15$ which is less than $1%$ higher than the Table I value and within the $7.5%$ error in Cf generated when the optimal parameter values are increased or decreased by one standard deviation. Given the relatively limited data near $βc=0$, the best approximation to the friction is generated by Eqs. (45) and (47). More data for weak pressure gradient conditions is needed.

In the presence of a pressure gradient, the wall parameters $(k,a,m)$ will be kept constant at the ZPG boundary layer values in Table I while the wake parameters $(b,n)$ will be treated as functions of βc using the correlations (45) and (47). This will be the approach used for the drag calculations presented in Secs. VI–VIII.

The important conclusion of this is that the length scale function (21) no longer depends on the parameters b and n independently, but depends only on the pressure gradient parameter βc.

(48)

and

(49)

It should be noted that the number of significant figures provided in this paper as well as in previous work on the UVP10,11,17 is intended to allow an interested reader to reproduce our results with accurate comparisons. It is not intended to reflect experimental or computed accuracy of the data.

In Sec. VI, the J0012 airfoil boundary layer calculation will be considered in detail but first it is necessary to analyze the flow at the airfoil leading edge. To model the leading edge flow depicted in Fig. 2, we will use the inviscid solution about an elliptic wing presented in Van Dyke.39 Let

(50)

where r is the leading edge radius of curvature, and c is the chord of the elliptic wing. The coordinate of the upper surface of the ellipse with the leading edge located at $x̃/c=0$ is

(51)

and the surface velocity is

(52)

The coordinate $x(x̃)$ along the wing surface is

(53)

where the function $E(θ,1−α)$ in Eq. (53) is the incomplete elliptic integral of the second kind.

Using Eqs. (51)–(53), the following limits are found to hold near the leading edge. Recall $ξ=x/r$

(54)

With the local inviscid solution (54) known, the viscous laminar solution near the leading edge can be generated from the classical Hiemenz solution at a stagnation point,

(55)

Substitute the variables in (55) into the vorticity equation

(56)

The result is

(57)

The fourth-order ordinary differential equation (ODE) in (57) can be integrated once leading to the third-order equation governing the Hiemenz solution.

(58)

The solution of (58) is shown in Fig. 20. Near the leading edge, the Hiemenz solution leads to the following limits:

(59)

where $f″(0)=1.232 59, c1=0.647 836$, and $c2=0.292 282$. The shape factor of the Hiemenz flow at the forward stagnation point is $H=δ1/δ2=2.216$, which can be compared to $H=2.5$ for channel flow and $H=2.604$ for the ZPG Blasius boundary layer.

FIG. 20.

Hiemenz streamwise velocity and vorticity.

FIG. 20.

Hiemenz streamwise velocity and vorticity.

Close modal

In the UVP integral boundary layer method presented in this paper, the simplest approach will be used where the UVP is assumed to apply all the way from the forward stagnation point of the airfoil. Therefore, it is important to consider the degree of error associated with this assumption by comparing the UVP with the “exact” Hiemenz solution.

The UVP approaches the laminar profile, (20), at the wing leading edge. From this profile, the following limits can be determined:

(60)

Substitute the limits in (60) into the Kármán equation (41). The result is

(61)

Equation (61) is a first-order ODE with the following integral:

(62)

where the initial condition $U(0)=0$ has been used. Substitute the last relation in (54), $U=αξ$, into (62). The result is

(63)

Compare the UVP limit (63) to the Hiemenz limit in Eq. (59). In the last relation in (59), let $ηh=2.381$ corresponding to $u/ue=0.99$,

(64)

The laminar limiting profile of the UVP, Eq. (20), can be expressed in terms of η by noting that $y/δh=η/ηh$

(65)

The Hiemenz and UVP velocity profiles are compared in Fig. 21.

FIG. 21.

Comparison between the Hiemenz velocity profile (black) and UVP (blue) near the forward stagnation point.

FIG. 21.

Comparison between the Hiemenz velocity profile (black) and UVP (blue) near the forward stagnation point.

Close modal

The conclusion from this discussion is that near the leading edge, the UVP limit (63) has the same dependence on α, Re, and ξ as the “exact” Hiemenz limit in (59) but the magnitude of $Rτ$, and therefore, the friction is about 28% below the Hiemenz limit. Keeping this in mind, the UVP will be assumed to apply beginning at the leading edge in the airfoil examples discussed later.

This is not the only approach that could be used. The laminar boundary layer equations could be solved beginning at the leading edge up to a transition region modeled with an appropriate transition criterion. This could then be followed by a turbulent boundary layer approximated by the UVP. For airfoil flows with a substantial laminar region, this would be more accurate but would require more information on free stream turbulence, and other factors that could influence transition. The viscous drag coefficient of low Reynolds number airfoil flows would be roughly $10%$ lower using this approach compared to using the UVP over the whole airfoil surface. The reason is that, even though the UVP friction might be lower very close to the leading edge as just discussed, the friction associated with the UVP channel flow limit is $10%$ or more higher than the Blasius-like friction of the rest of the laminar part of the flow,10 depending on the pressure gradient.

In this paper our focus is on high Reynolds number and so we will use the UVP all the way to the leading edge and comparisons will be with tripped airfoil data.

The J0012 airfoil geometry and inviscid surface velocity is shown in Fig. 22. Figure 23 shows the relation between the surface and chordwise coordinates including an inset showing the leading edge. The infinite slope of the wing near the leading edge and the singularity in the inviscid velocity derivative are both alleviated when the flow is expressed in terms of the surface coordinate $ξ=x/r$ rather than the chordwise coordinate $ξ̃=x̃/r$. The potential flow distributions of $U(ξ̃)$ for the J0012 and NACA0012 airfoils are shown in Fig. 24.

FIG. 22.

Joukowsky 0012 airfoil with free stream velocity and velocity derivative. Maximum surface velocity ($U=1.158 9$) is at $x̃/c=0.174 45$. Zero airfoil slope is at $x̃/c=0.368 47$.

FIG. 22.

Joukowsky 0012 airfoil with free stream velocity and velocity derivative. Maximum surface velocity ($U=1.158 9$) is at $x̃/c=0.174 45$. Zero airfoil slope is at $x̃/c=0.368 47$.

Close modal
FIG. 23.

Joukowsky 0012 surface coordinate vs chordwise coordinate. Inset shows shows the relation near the leading edge. Leading edge radius of curvature is $r/c=0.010 944$. Airfoil trailing edge is at $xc̃/r=91.3737$ corresponding to $xc/r=92.9941$.

FIG. 23.

Joukowsky 0012 surface coordinate vs chordwise coordinate. Inset shows shows the relation near the leading edge. Leading edge radius of curvature is $r/c=0.010 944$. Airfoil trailing edge is at $xc̃/r=91.3737$ corresponding to $xc/r=92.9941$.

Close modal
FIG. 24.

Potential flow velocity distributions about the J0012 (blue) and NACA0012 (magenta) airfoils.

FIG. 24.

Potential flow velocity distributions about the J0012 (blue) and NACA0012 (magenta) airfoils.

Close modal

Figure 25 shows eight successive solutions of the Kármán equation (41) over the J0012 airfoil at an airfoil chord Reynolds number, $Rchord=107$. The integral form of the UVP, Eq. (19), is used in this example. In Sec. VII, an alternative approach using the high Reynolds number explicit form of the UVP (23) will be described. The calculations were carried out in MathematicaTM on a 2017 MacBook Pro (OS 10.12.6, 2.9 GHz Intel core i7 processor). The iterative procedure converges relatively rapidly.

FIG. 25.

Successive iterations of $Rτ(x/r)$ for the J0012 airfoil. Inset shows the covergence of the trailing edge value of $RτTE$ with iteration number. The final iteration (blue) has a trailing edge value of $RτTE=4673$.

FIG. 25.

Successive iterations of $Rτ(x/r)$ for the J0012 airfoil. Inset shows the covergence of the trailing edge value of $RτTE$ with iteration number. The final iteration (blue) has a trailing edge value of $RτTE=4673$.

Close modal

Initial data required by the method begins with the kinematic viscosity ν and free stream velocity, $u∞$. The wing profile is defined, particularly the chord c and the leading edge radius r, which are used to determine α in Eqs. (50) and (63). The Reynolds number is $Re=u∞r/ν$. The wing profile is specified as

(66)

From the wing profile is calculated the wing surface coordinate x used as the independent variable in the UVP integral method,

(67)

See Figs. 23 and 24, where x and $x̃$ are normalized by the wing leading edge radius by multiplication by c/r. The velocity distribution $U(ξ)$, along with $dU/dξ$ are determined from a potential solution about the body in question. For the UVP integral method to work effectively, it is important that the free stream velocity and its derivative are accurately known and smooth.

Prior to integrating the Kármán equation or focusing on the airfoil, the integral form of the UVP, Eq. (19), is used to determine the functions $F0(Rτ), F1(Rτ), F2(Rτ)$, and $F3(Rτ)$ over the $Rτ$ range of interest. For the calculations presented here $0≤Rτ<1010$. This is done with the constants $(k,a,b,m,n)$ fixed at the values given in Table I, i.e., for $βc=0$. Since nested integration is required to produce the boundary layer thicknesses, this step can be quite slow when the range of interest includes large values of $Rτ$. However, these constant-b-constant-n thicknesses only need to be determined once and stored for use in the first iteration each time the drag for a new Rchord is to be computed.

In a bootstrap procedure, the Kármán equation (68) is solved numerically for the initial distribution $Rτ1(ξ)$ using the MathematicaTM function NDSolve,

(68)

The initial value of ξ is chosen to be close to the airfoil leading edge. For the J0012 case, with $U(ξ)$ known analytically, ξinitial was chosen to be $10−7$. For the NACA 0012 case discussed later, with $U(ξ)$ known approximately, ξinitial was chosen to be 0.1. The initial distribution, $Rτ1(ξ)$, is the black curve in Fig. 25 with $Rτ1TE=9784$.

The thickness functions $F0(Rτ1), F1(Rτ1), F2(Rτ1)$, and $F3(Rτ1)$ are computed over the range $0 and are used to calculate the first nonzero distribution of the modified Clauser parameter, $βc(Rτ1)$. After the first iteration, the wall parameters $(k,a,m)$ are kept constant at the ZPG boundary layer values in Table I while $n(βc)$ and $b(βc)$ are allowed to vary over the airfoil surface according to the correlations (45) and (47).

The monotonically increasing friction Reynolds number distribution over the airfoil computed in the first iteration, $Rτ1(ξ)$ together with $βc(Rτ1)$ are then used to prepare for iteration 2. New thickness functions $F01(Rτ1), F11(Rτ1), F21(Rτ1)$, and $F31(Rτ1)$ are computed over the range $0. The Kármán equation (69)

(69)

is integrated to generate $Rτ2(ξ)$. This is the red curve in Fig. 25.

The process is repeated until the nth iterate when further changes in $RτnTE$ are sufficiently small ($<2%$). This is the blue curve in Fig. 25. The iterative procedure is shown schematically in  Appendix B.

The time consuming steps in this procedure are the computations of the thickness functions, F1 and F2 which require nested integrations at npoints over the airfoil surface where npoints may vary from 500 for $Rchord=105$ up to 2500 for $Rchord=1012$. At $Rchord=1012$ each iteration may take 3–4 h on the laptop described above. In this example, eight iterations are needed to reach convergence with $Rτ1TE$ changing by less than $2%$ over the last two iterations. The corresponding change in $Cdv$ over the last two iterations is $0.3%$. The inset shows the convergence of $RτiTE(cr)$ at the end of which, $Rτ8TE=4673$. At this point, all of the boundary layer characteristics are known along with the velocity profile at any point.

The procedure could be repeated with the airfoil surface modified by the displacement thickness however the close correspondence between the NACA0012 potential flow and the boundary layer edge velocity from the SU2 converged solution at $Rchord=107$ shown in Fig. 35 indicates that further iterations would produce little change in the UVP drag.

Notice that the trailing edge values of βc in Fig. 29 never exceed $βc=18$. However, the user needs to be aware that during each iteration, the $βci$ calculated from the current $Rτi(ξ)$ distribution becomes very large, far larger than 18, as the trailing edge is approached where $uτ$ becomes relatively small (but not 0 because of the cusped trailing edge) and the displacement and momentum thickness Reynolds numbers in Eq. (43) become large. However note the decreasing values of$RτTE$with each iteration. It turns out that when $βci$ is used to calculate $Rτi+1$, the calculation reaches the trailing edge at a value of $Rτi+1TE$ that is considerably less than $RτiTE$ with a corresponding $βci$ that is less than 18.

Figures 26–32 show the convergence of the various functions required to integrate the Kármán equation (41) and determine the airfoil drag. Two other boundary layer properties of importance are the channel half height thickness (overall boundary layer thickness), δh, and the shape factor, $H=δ1/δ2$, shown in Figs. 32 and 33. The overall thickness is determined from

(70)
FIG. 26.

Iterative distributions of the wall friction coefficient Cf for the J0012 airfoil. Inset shows the covergence of the trailing edge value of $δh/c$ with iteration number. Final trailing edge value is $Cf=0.001 055$.

FIG. 26.

Iterative distributions of the wall friction coefficient Cf for the J0012 airfoil. Inset shows the covergence of the trailing edge value of $δh/c$ with iteration number. Final trailing edge value is $Cf=0.001 055$.

Close modal
FIG. 27.

Iterative distributions of the displacement thickness $δ1/c$ for the J0012 airfoil. Inset shows the covergence of the trailing edge value of $δ1/c$ with iteration number. Final trailing edge value is $δ1/c=0.006 201$.

FIG. 27.

Iterative distributions of the displacement thickness $δ1/c$ for the J0012 airfoil. Inset shows the covergence of the trailing edge value of $δ1/c$ with iteration number. Final trailing edge value is $δ1/c=0.006 201$.

Close modal
FIG. 28.

Iterative distributions of the momentum thickness $δ2/c$ for the J0012 airfoil. Inset shows the convergence of the trailing edge value of $δ2/c$ with iteration number. Final trailing edge value is $δ2/c=0.003 759$. All eight iterations overlay almost exactly.

FIG. 28.

Iterative distributions of the momentum thickness $δ2/c$ for the J0012 airfoil. Inset shows the convergence of the trailing edge value of $δ2/c$ with iteration number. Final trailing edge value is $δ2/c=0.003 759$. All eight iterations overlay almost exactly.

Close modal
FIG. 29.

Iterative distributions of the modified Clauser parameter, βc, for the J0012 airfoil. Inset shows the covergence of the trailing edge value of βc with iteration number. Final trailing edge value is $βc=17.238$. Note that $βc=0$ over the first iteration.

FIG. 29.

Iterative distributions of the modified Clauser parameter, βc, for the J0012 airfoil. Inset shows the covergence of the trailing edge value of βc with iteration number. Final trailing edge value is $βc=17.238$. Note that $βc=0$ over the first iteration.

Close modal
FIG. 30.

Iterative distributions of the UVP parameter b for the J0012 airfoil. Inset shows the covergence of the trailing edge value of b with iteration number. Final trailing edge value is $bTE=0.041 56$.

FIG. 30.

Iterative distributions of the UVP parameter b for the J0012 airfoil. Inset shows the covergence of the trailing edge value of b with iteration number. Final trailing edge value is $bTE=0.041 56$.

Close modal
FIG. 31.

Iterative distributions of the UVP parameter n for the J0012 airfoil. Inset shows the covergence of the trailing edge value of n with iteration number. Final trailing edge value is $nTE=6.0994$.

FIG. 31.

Iterative distributions of the UVP parameter n for the J0012 airfoil. Inset shows the covergence of the trailing edge value of n with iteration number. Final trailing edge value is $nTE=6.0994$.

Close modal
FIG. 32.

Iterative distributions of boundary layer thickness $δh/c$ for the J0012 airfoil. Inset shows the covergence of the trailing edge value of $δh/c$ with iteration number. Final trailing edge value is $δh/c=0.024 98$.

FIG. 32.

Iterative distributions of boundary layer thickness $δh/c$ for the J0012 airfoil. Inset shows the covergence of the trailing edge value of $δh/c$ with iteration number. Final trailing edge value is $δh/c=0.024 98$.

Close modal
FIG. 33.

Iterative distributions of the boundary layer shape factor $H=δ1/δ2$ for the J0012 airfoil. Inset shows the convergence of the trailing edge value of H with iteration number. Final trailing edge value is $H=1.6496$.

FIG. 33.

Iterative distributions of the boundary layer shape factor $H=δ1/δ2$ for the J0012 airfoil. Inset shows the convergence of the trailing edge value of H with iteration number. Final trailing edge value is $H=1.6496$.

Close modal

With the friction coefficient known over the airfoil surface, the viscous drag coefficient is determined using Eq. (8). The convergence of the J0012 and NACA0012 viscous drag coefficients at $Rchord=107$ are shown in Fig. 34. Figure 35 shows the potential flow over the NACA0012 wing.

FIG. 34.

Convergence of the J0012 (blue open circles) and NACA0012 (magenta open circles) viscous drag coefficients at $Rchord=107$. Converged values are $CdvJ0012=0.006 30$ and $CdvNACA0012=0.006 39$.

FIG. 34.

Convergence of the J0012 (blue open circles) and NACA0012 (magenta open circles) viscous drag coefficients at $Rchord=107$. Converged values are $CdvJ0012=0.006 30$ and $CdvNACA0012=0.006 39$.

Close modal
FIG. 35.

Potential flow velocity distribution about the NACA0012 airfoil (magenta) compared with the inviscid, boundary layer edge velocity (black) from SU2 at $Rchord=107$.

FIG. 35.

Potential flow velocity distribution about the NACA0012 airfoil (magenta) compared with the inviscid, boundary layer edge velocity (black) from SU2 at $Rchord=107$.

Close modal

In the presence of a pressure gradient, the wall parameters $(k,a,m)$ are kept constant at the ZPG boundary layer values in Table I while the wake parameters $(b,n)$ are treated as functions of βc using the correlations (45) and (47). Figure 36 shows the shape function, $ϕ$ Eq. (23), for a range of values of βc. Changes in $ϕ$ with βc are smooth and monotonic. A polynomial of high order is a safe, reliable fit in this situation since there is no question of $y/δh$ falling outside of the range $0≤y/δh≤1$. The family of shape functions shown in Fig. 36 is accurately approximated by the eighth-order polynomial

(71)
FIG. 36.

The shape function, $ϕ$, for various βc in the range $−1.0<βc<18.0$ evaluated at $Rτ=106$. Successive curves are in increments of 0.5.

FIG. 36.

The shape function, $ϕ$, for various βc in the range $−1.0<βc<18.0$ evaluated at $Rτ=106$. Successive curves are in increments of 0.5.

Close modal

The coefficients $c0(βc)$ to $c8(βc)$ are approximated by tenth-order polynomial functions of the pressure gradient parameter βc and are strictly limited to the range $−1.0<βc<18.0$ indicated in Fig. 36. The expressions for c0 to c8 are provided in  Appendix A. It is important to note that $ϕ(βc,y/δh)$ in Eq. (71) is universal in the sense that it applies to any pressure gradient wall flow that falls in the range $−1.0<βc<18.0$. Higher order fits could be used to increase the range of βc beyond 18.0.

The explicit high Reynolds number, variable pressure gradient form of the UVP reduces to the simple form

(72)

At high Reynolds number the parameter k, essentially the Kármán constant, is simply a scale factor on y+, u+ and the damping length scale a. Since b and n are correlated with the modified Clauser parameter, βc, the remaining optimal parameters governing the k-scaled UVP are only ka and m.

With $ϕ(ka,m,βc,y/δh)$ known explicitly as a bivariate polynomial, the integrals (26) and (27) can be carried out. This leads to polynomial expressions for the friction and thickness functions $F0(βc,Rτ), F1(βc,Rτ), F2(βc,Rτ)$, and $F3(βc,Rτ)$. Once these functions are known, the UVP integral method no longer requires the computation of nested integrals. As a result, the calculation time is independent of the Reynolds number. The viscous drag coefficient for any chord Reynolds number can be determined in a few seconds. Above about $Rchord=106$, this high Reynolds number explicit UVP method matches the results using the integral form of the UVP almost exactly. See Fig. 38 and Table IV.

TABLE IV.

Chord Reynolds number, $Cdv$ for the J0012 airfoil, $Cdv$ for the NACA0012 airfoil using both the integral and explicit forms of the UVP, $Cdv$ and $Cdp$ computed using SU2,33 tripped $Cdv$ data from Ladson.40

Rchord$CdvJ0012$$CdvNACA0012UVPintegral$$CdvNACA0012UVPexplicit$$CdvNACA0012SU2$$CdpNACA0012SU2$$CdvLadson80−180 grit$
105 0.014 797 0 0.014 817 4 0.013 499 8 0.012 746 0 0.004 738 7 ⋯
$5×105$ 0.010 352 9 0.010 397 7 0.010 032 5 0.010 108 7 0.002 664 7 ⋯
106 0.009 062 1 0.009 147 5 0.009 003 4 0.008 930 0 0.002 252 9 ⋯
$2×106$ 0.008 040 3 0.008 147 7 0.008 057 8 ⋯ ⋯ 0.008 53
$4×106$ 0.007 214 1 0.007 295 5 0.007 244 0 ⋯ ⋯ 0.007 33
$5×106$ 0.006 977 5 0.007 050 9 0.007 005 5 0.007 039 7 0.001 608 7 ⋯
$6×106$ 0.006 787 1 0.006 862 6 0.006 820 3 ⋯ ⋯ 0.006 82
$8.95×106$ 0.006 397 9 0.006 488 3 0.006 435 3 ⋯ ⋯ 0.006 51
107 0.006 296 5 0.006 394 3 0.006 294 3 0.006 431 9 0.001 495 7 ⋯
$1.2×107$ 0.006 137 3 0.006 228 17 0.006 170 8 ⋯ ⋯ 0.006 53
$5×107$ 0.005 112 8 0.005 102 1 0.005 082 9 0.005 128 0 0.001 092 0 ⋯
108 0.004 716 9 0.004 716 8 0.004 6502 0.004 696 2 0.000 933 7 ⋯
109 0.003 512 5 0.003 547 7 0.003 535 7 0.003 400 6 0.000 560 9 ⋯
1010 0.002 755 5 0.002 814 7 0.002 767 3 ⋯ ⋯ ⋯
1011 0.002 204 9 0.002 147 2 0.002 217 3 ⋯ ⋯ ⋯
1012 0.001 802 8 0.001 764 5 0.001 812 6 ⋯ ⋯ ⋯
Rchord$CdvJ0012$$CdvNACA0012UVPintegral$$CdvNACA0012UVPexplicit$$CdvNACA0012SU2$$CdpNACA0012SU2$$CdvLadson80−180 grit$
105 0.014 797 0 0.014 817 4 0.013 499 8 0.012 746 0 0.004 738 7 ⋯
$5×105$ 0.010 352 9 0.010 397 7 0.010 032 5 0.010 108 7 0.002 664 7 ⋯
106 0.009 062 1 0.009 147 5 0.009 003 4 0.008 930 0 0.002 252 9 ⋯
$2×106$ 0.008 040 3 0.008 147 7 0.008 057 8 ⋯ ⋯ 0.008 53
$4×106$ 0.007 214 1 0.007 295 5 0.007 244 0 ⋯ ⋯ 0.007 33
$5×106$ 0.006 977 5 0.007 050 9 0.007 005 5 0.007 039 7 0.001 608 7 ⋯
$6×106$ 0.006 787 1 0.006 862 6 0.006 820 3 ⋯ ⋯ 0.006 82
$8.95×106$ 0.006 397 9 0.006 488 3 0.006 435 3 ⋯ ⋯ 0.006 51
107 0.006 296 5 0.006 394 3 0.006 294 3 0.006 431 9 0.001 495 7 ⋯
$1.2×107$ 0.006 137 3 0.006 228 17 0.006 170 8 ⋯ ⋯ 0.006 53
$5×107$ 0.005 112 8 0.005 102 1 0.005 082 9 0.005 128 0 0.001 092 0 ⋯
108 0.004 716 9 0.004 716 8 0.004 6502 0.004 696 2 0.000 933 7 ⋯
109 0.003 512 5 0.003 547 7 0.003 535 7 0.003 400 6 0.000 560 9 ⋯
1010 0.002 755 5 0.002 814 7 0.002 767 3 ⋯ ⋯ ⋯
1011 0.002 204 9 0.002 147 2 0.002 217 3 ⋯ ⋯ ⋯
1012 0.001 802 8 0.001 764 5 0.001 812 6 ⋯ ⋯ ⋯

Referring back to Sec. II A, it should be pointed out that when Eq. (72) is used to generate the thicknesses, Eqs. (26) and (27), at high Reynolds number there is an error in neglecting the integration from the wall to the beginning of the log region. The log region for a ZPG boundary layer begins at approximately $y+=132$.17Figure 37 shows the ratio of the integration using the integral form of the UVP from the wall to $y+=132$ divided by integration using the explicit form of the UVP from 0 to $Rτ$, as a function of $Rτ$. Above $Rτ=105$ the error is less than $0.5%$.

FIG. 37.

Estimate of the error when the explicit form of the UVP, Eq. (72), is used to determine boundary layer thickness Reynolds numbers F1 (dotted line) and F2 (dashed line) neglecting the near wall region. Log region for a ZPG boundary layer begins at $y+=132$.17

FIG. 37.

Estimate of the error when the explicit form of the UVP, Eq. (72), is used to determine boundary layer thickness Reynolds numbers F1 (dotted line) and F2 (dashed line) neglecting the near wall region. Log region for a ZPG boundary layer begins at $y+=132$.17

Close modal

The J0012 airfoil was chosen for this study because it has an exact, and, therefore, highly accurate, potential flow solution and a well-defined cusped trailing edge. The NACA 0012 (with closed trailing edge) was chosen because this airfoil has been used as a standard for testing computational results, for assessing wall interference effects, and for developing correction procedures for wind tunnel testing. There is a variety of formulas in the literature for the NACA 0012 airfoil surface, consequently there is a variety of profiles used to study the flow, some with unclosed trailing edges and some with closed trailing edges. The profile used in this study is defined by

(73)

with constants, $a1=0.177 349 856, a2=−0.075 600 000, a3=−0.212 843 959 1, a4=0.173 640 303 0$, and $a5=−0.062 546 200 2$, which produce a closed trailing edge.

The NACA 0012 airfoil drag was determined using the UVP integral method and compared with computations of the drag using the aerodynamic design software suite, SU2.33 For the presentation, the boundary layer thicknesses, drag coefficient, and Reynolds number are normalized by the airfoil chord. Figure 35 shows the NACA0012 potential flow $U(ξ̃)$ determined using conformal mapping, compared with the free stream $U(ξ̃)$ computed from the SU2 pressure distribution at $Rchord=107$. The pressure drag from SU2 at this Reynolds number was $CdpSU2=0.001 496$.

The final Cdv vs Rchord results are shown in Fig. 38 and the data used to produce this figure are provided in Table IV. The nearly overlapping blue (J0012) and magenta (NACA0012) curves in Fig. 38 are generated using the UVP integral method with the mean wall parameter values, $(k¯,a¯,m¯)$, given in Table I. Recalling the discussion of Figs. 6 and 7 in Sec. II B, increasing or decreasing the optimal parameter values by one standard deviation would lead to a $7.5%$ change in Cf at $Rτ=5000$. The magenta error bar in Fig. 38 is approximately the variation of $Cdv$ that would occur if optimal parameters in the UVP were increased or decreased by the σ's given in Table I. The position of the magenta error bar at $Rchord=7×106$ is approximately at $Rτ=5000$. It should be noted that, while the error associated with changing the optimal parameters this way actually depends on $Rτ$, the magenta error bar provides a reasonable approximation to the error in the calculated $Cdv$ based on the set of boundary layer cases evaluated to date.

FIG. 38.

$Cdv$ for the J0012 airfoil (blue curve, see inset) and NACA 0012 airfoil (magenta curve) determined using the integral form of the UVP. The magenta vertical error bar is the variation of $Cdv$ that would occur if optimal parameters in the UVP were increased or decreased by the variances given in Table I. The green curve is $Cdv$ calculated using the high Reynolds number explicit form of the UVP. Red filled circles are $Cdv$ computed using SU233 with the Spalart–Allmaras model.37 Open black circles are NACA 0012 tripped data from Ladson40 averaged over grit sizes 80–180. The black dashed curve is the correlation of tripped NACA 0012 data from McCroskey.41 The black vertical error bar represents the variation ($ΔCdv=0.001 2$) in the tripped data evaluated by McCroskey.

FIG. 38.

$Cdv$ for the J0012 airfoil (blue curve, see inset) and NACA 0012 airfoil (magenta curve) determined using the integral form of the UVP. The magenta vertical error bar is the variation of $Cdv$ that would occur if optimal parameters in the UVP were increased or decreased by the variances given in Table I. The green curve is $Cdv$ calculated using the high Reynolds number explicit form of the UVP. Red filled circles are $Cdv$ computed using SU233 with the Spalart–Allmaras model.37 Open black circles are NACA 0012 tripped data from Ladson40 averaged over grit sizes 80–180. The black dashed curve is the correlation of tripped NACA 0012 data from McCroskey.41 The black vertical error bar represents the variation ($ΔCdv=0.001 2$) in the tripped data evaluated by McCroskey.

Close modal

The green (NACA0012) curve was generated using the explicit form of the UVP which allows high Reynolds number cases to be computed in a few seconds. The three UVP curves are almost identical at Reynolds numbers above $Rchord>106$.

The open black circles in Fig. 38 are from Ladson40 who used carborundum strips placed at $x̃/c=0.05$ to generate a tripped boundary layer over an NACA0012 wing. The data points shown are averages over grit sizes 80, 120, and 180 and Mach numbers 0.15 and 0.30. The black dashed curve in Fig. 38 is a fit to high quality tripped data from several sources (including Ladson) collected and evaluated by McCroskey41

(74)

Note that an interpolation of the pressure drag coefficient, $CdpSU2$ is subtracted from the Cd data in Ladson40 and McCroskey41 to generate the black open circles and dashed line in Fig. 38.

The McCroskey correlation (74) lies slightly below the UVP results and has a slightly smaller slope and intersects the blue, magenta, and green curves at about $Rchord=107$. The black error bar in Fig. 38 at $Rchord=107$ indicates the variation in the data in Fig. 4 of McCroskey.41

The filled red circles in Fig. 38 are results from SU233 using the Spalart–Allmaras model.37 No transition model was used in the computations. Between $Rchord=106$ and 109, the difference between the UVP drag and the SU2 data are generally less than 0.0002. The results for $Rchord=109$ were obtained using a procedure consistent with the lower Rchord cases. The finest grid available from the Langley Turbulence Modeling Resource webpage42 gave a wall spacing of $10−7$ m corresponding to a minimum $y+≳1$ over significant portions of the airfoil instead of the preferred $y+≤1$. This inability to accurately capture the viscous sublayer over much of the airfoil may explain the slight under-prediction of SU2 relative to the UVP at this high Re. While further validation in high Re flows is still needed, this case does demonstrate close agreement of the UVP method with established turbulence models at high Re.

The good agreement between the UVP and SU2 results is to be expected given the fact that the flat plate wall friction coefficients generated by both Eqs. (30) and (34) are in good agreement.

The results in Fig. 38 are all consistent with the point made in Sec. II C that the UVP represents a tripped boundary layer with a critical Reynolds number close to the minimum needed to sustain turbulent flow.

To investigate the UVP and SU2 drag coefficient results in a little more detail, the viscous drag force distribution along the NACA 0012 airfoil is shown in Figs. 39 and 40.

FIG. 39.

Viscous drag force on the NACA0012 airfoil at $Rc=107$ as a function of position. Comparison between the UVP first iteration (magenta, dashed) and the UVP final iteration (magenta, solid).

FIG. 39.

Viscous drag force on the NACA0012 airfoil at $Rc=107$ as a function of position. Comparison between the UVP first iteration (magenta, dashed) and the UVP final iteration (magenta, solid).

Close modal
FIG. 40.

Viscous drag force on the NACA0012 airfoil at $Rc=107$ as a function of position. Comparison between the UVP final iteration (magenta) and the viscous drag computed using SU2 with the Spalart–Allmaras model (black).

FIG. 40.

Viscous drag force on the NACA0012 airfoil at $Rc=107$ as a function of position. Comparison between the UVP final iteration (magenta) and the viscous drag computed using SU2 with the Spalart–Allmaras model (black).

Close modal

Figure 39 shows the drag force distribution generated by the first and last iterations of the UVP integral method. The first iteration is remarkably good considering that there is no adjustment of b or n. Most of the variation in drag is accounted for by the relationship between the free stream velocity and wall friction inherent to the UVP. The largest difference between the first and last iterations occurs near the airfoil trailing edge where the adverse pressure gradient adjustments of b and n are needed to properly capture the drag force as well as other boundary layer characteristics in that region.

In Fig. 40 the SU2 distribution, shown in black, is compared to the final iteration of the UVP distribution, shown in magenta. The two curves almost overlap over most of the wing. Over the last 2% of the chord the UVP drag drops rapidly to almost zero near the finite angle trailing edge, while the SU2 curve rises slightly in response to the beginning of the airfoil wake.

The universal velocity profile, developed originally for pipe flow,11 provides a very accurate approximation to channel flow and boundary layer flow with zero, favorable, and adverse pressure gradient. A modified Clauser parameter, βc, is defined for the first time and shown in Figs. 18 and 19 to correlate well with the parameters b and n that characterize the wake portion of the UVP. This result provides a complete characterization of the Reynolds number and pressure gradient dependence of the UVP enabling it to be used to create the new integral method presented in this paper. Interestingly, the first iteration on Cdv with b and n fixed at the zero pressure gradient values is remarkably close to the final value, although the distributions of $Rτ$ and Cf generated by the modified Clauser parameter are essential to properly capture the trailing edge flow.

For all calculations, the parameters that characterize the wall flow, $(k¯,a¯,m¯)$ were held constant at the values in Table I independent of the pressure gradient.17 As additional boundary layer data become available, the average optimal parameter values will continue to change slightly although the impact on calculated values of $Cdv$ will be decreasingly small as the sample size grows.

The $Cdv$ generated by the UVP in Fig. 38 and Table IV is generally in very good agreement with the McCroskey/Ladson data and the $Cdv$ data generated by SU2. The conclusion that can be drawn from the results in Table I and Fig. 38 is that the UVP integral method is an effective method for determining viscous drag on smooth aerodynamic shapes. The conclusion from the discussion of Figs. 5 and 9 and the comparison between the UVP and Ladson/McCroskey data in Fig. 38 is that when the UVP is used beginning at the airfoil leading edge, the drag should be regarded as that of a tripped boundary layer with a critical Reynolds number close to the minimum required to sustain turbulent flow.

To handle an untripped case, the laminar and transitional part of the flow will need to be determined using a transition criterion for the boundary layer along with additional information about the free stream. One might use the well-known model of boundary layer transition where transition is considered complete after the amplification of unstable initial disturbances by a factor eN.43,44 For a boundary layer developing in a very low free stream turbulence environment, for example, Schubauer and Klebanoff,30N would be selected to be in the range 9–12, whereas for a boundary layer developing in a highly turbulent free stream one might choose $N≤1$. The problem of how best to incorporate a transition criterion into the universal velocity profile is a subject of current research.

Large eddy simulation of wall flows is faced with the conflict between the need to reduce grid resolution near the wall in order to reduce computational cost and the requirement to accurately determine the wall friction. Since the large eddy simulation (LES) grid under-resolves the fine scale motions near the wall that are responsible for generating the wall stress, empirical wall functions need to be used between the first grid point and the wall in order to determine the friction. These functions inevitably encounter the log-layer mismatch problem.45–48 The UVP provides an accurate value of the mean stream-wise and wall-normal velocity components as well as the turbulent shear stress throughout the viscous wall layer without assuming a log or power law profile and may provide a useful alternative to the current approach. That the UVP can provide an accurate connection to the wall is evidenced by the fact that virtually all the experimental non-DNS datasets that we have approximated thus far only contain data that is well outside the wall layer.

The inherent dependence of the UVP on Reynolds number, extended to include the effect of pressure gradient enables it to be used as the basis of a new method for integrating the Kármán equation for a wide variety of attached, wall bounded flows. There is really no practical limit to the Reynolds number that can be computed suggesting that, with modification for the effects of compressibility and/or roughness, it can be applied to very large-scale aerodynamic, hydrodynamic, and geophysical flows.

Ultimately, fundamental questions about the optimal parameters will need to be answered as to their numerical values, dependence on flow geometry (seen clearly in Table I), and possible weak dependence on Reynolds number, free stream turbulence, surface roughness, Mach number, and so on. It is not clear when the high Reynolds number boundary layer data required to address these questions will become available.

The authors would like to express appreciation to Stanford University for continued support. The computational resources made available through the Stanford Sherlock High Performance Computing Cluster are gratefully acknowledged. All grids used for the computations came from the Langley Turbulence Modeling Resource webpage.

The authors have no conflicts to disclose.

Brian Joseph Cantwell: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (lead); Investigation (equal); Methodology (lead); Resources (equal); Software (equal); Supervision (lead); Writing – original draft (lead); Writing – review & editing (lead). Eylul Bilgin: Conceptualization (supporting); Data curation (supporting); Formal analysis (supporting); Investigation (equal); Methodology (supporting); Resources (equal). Jacob Troy Needels: Investigation (equal); Methodology (supporting); Resources (equal); Software (equal).

The data that support the findings of this study are available within the article.

The polynomial approximations to the coefficients in Eq. (71) are listed below. These approximations are only valid in the range $−1.0<βc<18.0$. Outside that range the higher order terms will produce very large errors,

(A1)
(A2)
(A3)
(A4)
(A5)
(A6)
(A7)
(A8)
(A9)

A flow chart showing the procedure for the UVP integral boundary layer method is shown in Fig. 41.

FIG. 41.

Flowchart summarizing the UVP method.

FIG. 41.

Flowchart summarizing the UVP method.

Close modal
1.
L.
Prandtl
, “
The mechanics of viscous fluids
,” in
Aerodynamic Theory
, edited by
W. F.
Durand
(
Guggenheim Fund for the Promotion of Aeronautics
,
1934
), Vol.
III
, pp.
80
90
.
2.
B.
Thwaites
,
Incompressible Aerodynamics
(
Oxford University Press
,
1960
), pp.
62
64
.
3.
K.
Pohlhausen
, “
Zur naherungsweisen integration der differentialgleichungen der laminaren reilbungsschicht
,”
ZAMM
1
,
252
260
(
1921
).
4.
J.
Majdalani
and
L.
Xuan
, “
On the Kármán momentum-integral approach and the Pohlhausen paradox
,”
Phys. Fluids
32
,
123605
(
2020
).
5.
M. R.
, “
Entrainment in the turbulent boundary layer
,” Aeronautical Research Council Reports and Memoranda No. 3152,
1960
.
6.
H.
Ludwieg
and
W.
Tillman
, “
Investigations of wall shearing stress in turbulent boundary layers
,” NACA Technical Memorandum No. 1285,
1950
.
7.
H.
Schlichting
and
K.
Gersten
,
Boundary Layer Theory
(
Springer
,
2000
).
8.
B. J.
Cantwell
,
Introduction to Symmetry Analysis
(
Cambridge University Press
,
2002
).
9.
T.
von Kármán
, “
Uber laminaire und turbulente reibung
,”
ZAMM
1
,
233
252
(
1921
).
10.
B. J.
Cantwell
, “
Integral measures of the zero pressure gradient boundary layer over the Reynolds number range $0≤Rτ<∞$
,”
Phys. Fluids
33
,
085108
(
2021
).
11.
B. J.
Cantwell
, “
A universal velocity profile for smooth wall pipe flow
,”
J. Fluid Mech.
878
,
834
(
2019
).
12.
T.
von Kármán
, “
Mechanical similitude and turbulence
,” NACA Technical Memorandum No. 611,
1931
.
13.
L.
Prandtl
, “
The mechanics of viscous fluids
,” in
Aerodynamic Theory
, edited by
W. F.
Durand
(
Guggenheim Fund for the Promotion of Aeronautics
,
1934
), Vol III, pp.
102
112
.
14.
E. R.
van Driest
, “
On turbulent flow near a wall
,”
J. Aeronaut. Sci.
23
,
1007
1011
(
1956
).
15.
W.
Kollmann
, “
Asymptotic properties of mixing length closures for turbulent pipe flow
,”
Phys. Fluids
32
,
115126
(
2020
).
16.
T.
Cebeci
and
P.
,
Momentum Transfer in Boundary Layers
(
Hemisphere Publishing Corporation
,
1977
), Chap. 6.4, p.
174
.
17.
M. A.
Subrahmanyam
,
B. J.
Cantwell
, and
J. J.
Alonso
, “
A universal velocity profile for turbulent wall flows including adverse pressure gradient boundary layers
,”
J. Fluid Mech.
933
,
A16
(
2022
).
18.
D. I.
Pullin
,
M.
Inoue
, and
N.
Saito
, “
On the asymptotic state of high Reynolds number, smooth-wall turbulent flows
,”
Phys. Fluids
25
,
015116
(
2013
).
19.
M. A.
Subrahmanyam
,
B. J.
Cantwell
, and
J. J.
Alonso
, “
A universal velocity profile for turbulent wall flows
,” AIAA Paper No. 2021-0061,
2021
.
20.
M. P.
Simens
,
J.
Jiménez
,
S.
Hoyas
, and
Y.
Mizuno
, “
A high-resolution code for turbulent boundary layers
,”
J. Comput. Phys.
228
,
4218
4231
(
2009
).
21.
G.
Borrell
,
J. A.
Sillero
, and
J.
Jiménez
, “
A code for direct numerical simulation of turbulent boundary layers at high Reynolds numbers in BG/P supercomputers
,” in
23rd International Conference on Parallel Fluid Dynamics ParCFD2011, 2011
[
Comput. Fluids
80
,
37
43
(
2013
)].
22.
J. A.
Sillero
,
J.
Jiménez
, and
R. D.
Moser
, “
One-point statistics for turbulent wall-bounded flows at Reynolds numbers up to $δ+$ = 2000
,”
Phys. Fluids
25
,
105102
(
2013
).
23.
G.
Eitel-Amor
,
O.
Ramis
, and
P.
Schlatter
, “
Simulation and validation of a spatially evolving turbulent boundary layer up to Re = 8300
,”
Int. J. Heat Fluid Flow
47
,
57
69
(
2014
).
24.
P. S.
Klebanoff
, “
Characteristics of turbulence in a boundary layer with zero pressure gradient
,”
NACA Report No. 1247
,
1955
.
25.
D. B.
DeGraaff
and
J. K.
Eaton
, “
Reynolds-number scaling of the flat-plate turbulent boundary layer
,”
J. Fluid Mech.
422
,
319
346
(
2000
).
26.
R.
Baidya
,
J.
Phillip
,
N.
Hutchins
,
J. P.
Monty
, and
I.
Marusic
, “
Distance from the wall scaling of turbulent motions in wall-bounded flows
,”
Phys. Fluids
29
,
020712
(
2017
).
27.
A. E.
Perry
and
I.
Marusic
, “
A wall-wake model for the turbulence structure of boundary layers. Part 2. Further experimental support
,”
J. Fluid Mech.
298
,
361
407
(
1995
).
28.
M.
Jones
,
I.
Marusic
, and
A. E.
Perry
, “
Evolution and structure of sink-flow turbulent boundary layers
,”
J. Fluid Mech.
428
,
1
27
(
2001
).
29.
M. V.
Zagarola
and
A. J.
Smits
, “
Mean-flow scaling of turbulent pipe flow
,”
J. Fluid Mech.
373
,
33
79
(
1998
).
30.
G. B.
Schubauer
and
P. S.
Klebanoff
, “
Contribution on the mechanics of boundary layer transition
,”
NACA Technical Note TN 3489
,
1955
.
31.
J.
Coupland
, “
ERCOFTAC special interest group on laminar to turbulent transition and retransition T3A and T3B test cases
,”
Technical Report No. A309514
,
1990
.
32.
U.
Kaynak
,
O.
Bas
,
S. C.
Cakmakcioglu
, and
I. H.
Tuncer
, “
Transition modeling for low to high speed boundary layer flows with CFD applications
,” in
Boundary Layer Flows—Theory, Applications, and Numerical Methods
(
IntechOpen
,
2020
).
33.
T. D.
Economon
,
F.
Palacios
,
S.
Copeland
,
T.
Lukaczyk
, and
J. J.
Alonso
, “
SU2: An open-source suite for multi-physics simulation and design
,”
AIAA J.
54
,
828
846
(
2016
).
34.
J. F.
Nash
, “
A note on skin-friction laws for the incompressible turbulent boundary layer
,” Aeronautical Research Council Current Papers No. 862,
1966
.
35.
H.
Nagib
,
C.
Christophorou
,
J.
Reudi
,
P.
Monkewitz
,
J.
Osterlund
, and
S.
Gravante
, “
Can we ever rely on results from wall-bounded turbulent flows without direct measurements of wall shear stress?
,” AIAA Paper No. 2004-2392,
2004
.
36.
M.
Polewski
and
P. G. A.
Cizmas
, “
Several cases for the validation of turbulence models implementation
,”
J. Appl. Sci.
11
,
3377
(
2021
).
37.
P. R.
Spalart
and
S. R.
Allmaras
, “
A one equation turbulence model for aerodynamic flows
,” AIAA Paper No. 89-0366,
1992
.
38.
F. H.
Clauser
, “
,”
J. Aeronaut. Sci.
21
,
91
108
(
1954
).
39.
M.
Van Dyke
,
Perturbation Methods in Fluid Mechanics
(
,
1964
).
40.
C. L.
, “
Effects of independent variation of Mach and Reynolds numbers on the low-speed aerodynamic characteristics of the NACA 0012 airfoil section
,” NASA Technical Memorandum No. 4074,
1988
.
41.
W. J.
McCroskey
, “
A critical assessment of wind tunnel results for the NACA 0012 airfoil
,” NASA Technical Memorandum No. 100019,
1987
.
42.
D. C.
Jesperson
,
T.
Pulliam
, and
M. L.
Childs
, “
OVERFLOW turbulence modeling resource validation results
,” “
NASA Technical Report No. NAS-2016-01
,
2016
.
43.
A. M. O.
Smith
and
N.
Gamberoni
, “
Transition, pressure gradient and stability theory
,”
Douglas Aircraft Report No. ES 26388
,
1956
.
44.
J. L.
van Ingen
, “
A suggested semi-empirical method for the calculation of the boundary layer transition region
,”
Technische Hogeschool Vliegtuigbouwkunde Report No. VTH-74
,
1956
.
45.
J.
Brasseur
and
T.
Wei
, “
Designing large-eddy simulation of the turbulent boundary layer to capture law-of-the-wall scaling
,”
Phys. Fluids
22
,
021303
(
2010
).
46.
X.
Yang
,
G.
Park
, and
P.
Moin
, “
Log-layer mismatch and modeling of the fluctuating wall stress in wall-modeled large-eddy simulations
,”
Phys. Rev. Fluids
2
,
104601
(
2017
).
47.
S.
Kawai
and
J.
, “
Wall-modeling in large eddy simulation: Length scales, grid resolution, and accuracy
,”
Phys. Fluids
24
,
015105
(
2012
).
48.
S.
Cai
and
P.
Sagaut
, “
Explicit wall models for large eddy simulation
,”
Phys. Fluids
33
,
041703
(
2021
).