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.

## I. INTRODUCTION

In volume III of the Durand series on Aerodynamic Theory, Prandtl^{1} 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.

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 Thwaites^{2} 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)

subject to the no-slip and free stream conditions

where *u*, *u _{e}*, 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).

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\u2032v\u2032\xaf$, the calculation of boundary layer characteristics can proceed.

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

The friction velocity is defined as

The displacement thickness is defined as

and the momentum thickness is

where *δ _{h}* indicated in Fig. 1 is a boundary layer thickness defined in Cantwell

^{10}as the equivalent channel half height thickness. The thickness,

*δ*, will be discussed further in Sec. II B.

_{h}Wing drag is calculated using the construction in Fig. 3. The differential force acting in the $x\u0303$ direction on, say, the wing upper surface is

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

where the local (based on *u _{e}*) friction coefficient is

## II. THE UNIVERSAL VELOCITY PROFILE (UVP)

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

where the pressure gradient is a constant, Cantwell.^{11}

Express Eq. (10) in wall normalized coordinates

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

where $R\tau $ is the friction Reynolds number

The turbulent shear stress is modeled using classical mixing length theory (von Kármán;^{12} see also Prandtl^{13} and van Driest^{14}). Let

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

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

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

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

independent of the choice of $\lambda (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 $\lambda (y+)$ at a given $R\tau $

where $u+(0)=0$ and $u+(R\tau )=ue/u\tau $. At low Reynolds number, $\lambda \u21920$ and Eq. (19) approaches the laminar pipe/channel flow solution

where, in the laminar limit, $R\tau =(2ue\delta h/\nu )1/2$.

The mixing length model introduced by Cantwell^{11} to approximate pipe data is

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 Driest^{14} except for the exponent *m* that helps to determine the shape and thickness of the near wall velocity profile. Near the wall

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\tau >2000/k$),^{11} *k* 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\u2264y\u2264\delta h$ and $0\u2264R\tau <\u221e$. 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.

### A. The high Reynolds number explicit form of the UVP

The classical wall-wake formulation is silent as to the Reynolds number at which the formulation actually applies. Whereas for $R\tau \u22652000/k$, the integral form of the UVP, Eq. (19) approaches the explicit form^{10,11}

Importantly, the shape function, $\varphi $, is independent of $R\tau $. As long as $R\tau $ 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 $\varphi $. 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\tau \u22652000/k$, say $R\tau =105$ or 10^{6}, to produce $\varphi $. Any velocity profile with $R\tau \u22652000/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\tau $ is generated using Eq. (23). Complete details are given in Cantwell^{10,11} and Subrahmanyam *et al.*^{17}

Flow . | $k\xaf$ . | σ
. _{k} | $a\xaf$ . | σ
. _{a} | $m\xaf$ . | σ
. _{m} | $b\xaf$ . | σ
. _{b} | $\u2009\u2009\u2009n\xaf$ . | σ
. _{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\xaf$ . | σ
. _{k} | $a\xaf$ . | σ
. _{a} | $m\xaf$ . | σ
. _{m} | $b\xaf$ . | σ
. _{b} | $\u2009\u2009\u2009n\xaf$ . | σ
. _{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 |

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

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\tau $. Note that the two curves merge around $R\tau =2000/k$.

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 $\u2202u/\u2202y=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\u2264R\tau <\u221e$. 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; *F*_{0}, *F*_{1}, *F*_{2}, and *F*_{3}. As noted above, the boundary layer friction law is generated by evaluating the UVP at $y+=R\tau $

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

and the local momentum thickness Reynolds number is written as

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

### B. Optimal parameters

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

Parameter values $(k\xaf,a\xaf,m\xaf,b\xaf,n\xaf)$ 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\tau $ 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\xaf,a\xaf,m\xaf,b\xaf,n\xaf)=(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 data^{11} 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\tau $ 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\tau =5000$, the friction coefficient based on the mean parameters in Table I is $Cf=0.002\u2009378$. With the mean parameters all increased by one standard deviation $Cf+=0.002\u2009463$ and with the parameters decreased by one standard deviation $Cf\u2212=0.002\u2009277$ for an overall variation of approximately $7.5%$. In Sec. VIII where airfoil drag comparisons are discussed, this uncertainty in *C _{f}* will be discussed again.

x (m)
. | $ue(ms)$ . | $R\delta 1$ . | $R\delta 2$ . | β
. | β
. _{c} | $\delta 998(m)$ . | $R\delta 998$ . | $R\tau $ . | $(ueu\tau )$ . | k
. | a
. | m
. | b
. | n
. | $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\delta 1$ . | $R\delta 2$ . | β
. | β
. _{c} | $\delta 998(m)$ . | $R\delta 998$ . | $R\tau $ . | $(ueu\tau )$ . | k
. | a
. | m
. | b
. | n
. | $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 |

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 Cantwell^{10} 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.

#### 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

*δ*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.

_{h}^{10,11,17}

The UVP has a well-defined outer edge, where $u/ue=1$ and $\u2202u/\u2202y=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 $\u2202u/\u2202y$ is slightly greater than zero. Probably the most common choice is the boundary layer thickness corresponding to $ue/Ufreestream=0.99$, where *U _{freestream}* 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

*δ*is too small, the data will be cutoff too far short of the boundary layer edge. The condition $\u2202u/\u2202y=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 Cantwell

_{h}^{10}where

*δ*is first defined.

_{h}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\u22120.998)Ufreestream$.

### C. Evolution of the UVP from the laminar to the turbulent profile as $R\tau $ increases from 0 to $\u221e$

Figure 8 shows how the shape of the universal velocity profile becomes fuller as $R\tau $ 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}

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\tau =30$ corresponding to $R\delta 1=122,\u2009R\delta 2=51$ and $Rx=$15 700 with a shape factor *H* = 2.39. The profile approaches a fully turbulent shape at about $R\tau =500$ corresponding to $R\delta 1=2030,R\delta 2=1373$ and $Rx=645\u2009000$ 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\tau =121,\u2009R\delta 1=520$, and $Rx=104\u2009000$,^{7} to a nominal transition Reynolds number, $R\delta 1=1600,\u2009R\tau =392$ and $Rx=474\u2009000$. 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 *T _{u}*. 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 Klebanoff

^{30}and Coupland.

^{31}Note that large $ue/u\tau $ implies low

*C*. Extremely low values of

_{f}*T*can lead to critical Reynolds numbers exceeding $R\tau =2000$ corresponding to approximately $Rex=2\u2009700\u2009000$ while the UVP implies turbulent values of the friction coefficient beginning at about $R\tau \u2248200$ corresponding to $Rex\u2248200\u2009000$.

_{u}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 *T _{u}*. 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\tau transition\u2248200$. In other words, the UVP corresponds to a tripped boundary layer.

Figure 9 can be compared to Fig. 14 in Cantwell^{11} 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.

### D. Friction law comparisons

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\tau =2000/k$, the UVP friction law is

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

where *H* is the shape factor, $H=\delta 1/\delta 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:

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

Finally, Polewski and Cizmas^{36} 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–Allmaras^{37} 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

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

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<R\tau <10\u2009000$, the Ludwieg–Tillman law leads to unrealistically high $ue/u\tau $ (low *C _{f}*) 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<R\tau <25\u2009000$ which covers the vast bulk of available boundary layer data ($0<R\delta 1<1.009\xd7105$ corresponding to $0<R\delta 2<7.93\xd7104$ corresponding to $0<Rx<7.38\xd7107$). All five curves in Fig. 10 agree within about 5% over the range above $R\tau >1000$. The closest agreement with the UVP law in this Reynolds number range are the laws due to Nash^{34} and Spalart and Allmaras.^{37} Note that only the UVP merges with the laminar limit, $R\tau \u21920$ 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\u2009000<R\delta 2<70\u2009000$. 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 Nash^{34} (magenta dashed line). Note that the $R\tau $ and $R\delta 2$ domains in Figs. 11 and 12 represent almost identical Reynolds number ranges.

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

## III. THE VON KÁRMÁN INTEGRAL EQUATION EXPRESSED IN TERMS OF $R\tau $

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

The main drag results are expressed in terms of the chord Reynolds number, *R _{chord}*, 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 $\xi =x/r$. The dimensionless potential flow velocity at the surface of the body is

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

Differentiate (37) with respect to *x*

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.

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

## IV. THE EFFECT OF PRESSURE GRADIENT ON THE UVP

Subrahmanyam *et al.*^{17} used the UVP to approximate the adverse pressure gradient data of Perry and Marusic^{27} for the two upstream velocity cases, $u\u221e=10$ and $u\u221e=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.

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,

where $u\u221e$ 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\tau $ that characterize the data. Comparable values of *k* were observed at low $R\tau $ in the pipe data (Cantwell^{11} Fig. 9, PSP cases 1–5). The mean parameter values in Table I for pipe flow were averaged over PSP cases 6–26.

$K\xd7107$ . | x (mm)
. | $u0(ms)$ . | $R\delta 1$ . | $R\delta 2$ . | $\u2212\beta $ . | $\u2212\beta c$ . | $R\delta 996$ . | $R\tau $ . | $(ueu\tau )$ . | k
. | a
. | m
. | b
. | n
. | $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\xd7107$ . | x (mm)
. | $u0(ms)$ . | $R\delta 1$ . | $R\delta 2$ . | $\u2212\beta $ . | $\u2212\beta c$ . | $R\delta 996$ . | $R\tau $ . | $(ueu\tau )$ . | k
. | a
. | m
. | b
. | n
. | $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 |

### A. A modified Clauser parameter

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 Clauser^{38} parameter $\beta =(\delta 1/\tau w)(dp/dx)$. See, for example, Perry and Marusic.^{27} Here, we will introduce a modified Clauser parameter

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)

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

with rms error 0.0440. Equation (46) gives a good approximation to most of the data but is not accurate near $\beta 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

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

The effect of this change can be seen by evaluating *C _{f}* at $\beta c=0$ for $R\tau =104$. If the parameters

*b*and

*n*from Table I, $(b,n)=(0.1752,2.1707)$, are used $Cf=0.002\u200913$. If Eqs. (45) and (46) are used, $(b,n)=(0.3050,1.4194)$ and $Cf=0.002\u200938$, 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\u200915$ which is less than $1%$ higher than the Table I value and within the $7.5%$ error in

*C*generated when the optimal parameter values are increased or decreased by one standard deviation. Given the relatively limited data near $\beta c=0$, the best approximation to the friction is generated by Eqs. (45) and (47). More data for weak pressure gradient conditions is needed.

_{f}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}*.

and

It should be noted that the number of significant figures provided in this paper as well as in previous work on the UVP^{10,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.

## V. LEADING EDGE TREATMENT

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

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\u0303/c=0$ is

and the surface velocity is

The coordinate $x(x\u0303)$ along the wing surface is

where the function $E(\theta ,1\u2212\alpha )$ 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 $\xi =x/r$

### A. The Hiemenz solution near the leading edge

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,

Substitute the variables in (55) into the vorticity equation

The result is

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

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

where $f\u2033(0)=1.232\u200959,\u2009c1=0.647\u2009836$, and $c2=0.292\u2009282$. The shape factor of the Hiemenz flow at the forward stagnation point is $H=\delta 1/\delta 2=2.216$, which can be compared to $H=2.5$ for channel flow and $H=2.604$ for the ZPG Blasius boundary layer.

### B. The UVP near the leading edge

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:

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

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

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

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

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

The conclusion from this discussion is that near the leading edge, the UVP limit (63) has the same dependence on *α*, *R _{e}*, and

*ξ*as the “exact” Hiemenz limit in (59) but the magnitude of $R\tau $, 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.

## VI. ITERATIVE PROCEDURE AND SOLUTION FOR THE J0012 AIRFOIL AT $Rchord=107$ USING THE INTEGRAL FORM OF THE UVP

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 $\xi =x/r$ rather than the chordwise coordinate $\xi \u0303=x\u0303/r$. The potential flow distributions of $U(\xi \u0303)$ for the J0012 and NACA0012 airfoils are shown in Fig. 24.

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 *Mathematica*^{TM} on a 2017 MacBook Pro (OS 10.12.6, 2.9 GHz Intel core i7 processor). The iterative procedure converges relatively rapidly.

### A. Step 1—Initial data, initial calculation of $R\delta 1(R\tau )$ and $R\delta 2(R\tau )$, first integration of the Kármán equation

Initial data required by the method begins with the kinematic viscosity *ν* and free stream velocity, $u\u221e$. 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\u221er/\nu $. The wing profile is specified as

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

See Figs. 23 and 24, where *x* and $x\u0303$ are normalized by the wing leading edge radius by multiplication by *c*/*r*. The velocity distribution $U(\xi )$, along with $dU/d\xi $ 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\tau ),\u2009F1(R\tau ),\u2009F2(R\tau )$, and $F3(R\tau )$ over the $R\tau $ range of interest. For the calculations presented here $0\u2264R\tau <1010$. This is done with the constants $(k,a,b,m,n)$ fixed at the values given in Table I, i.e., for $\beta 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\tau $. 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 *R _{chord}* is to be computed.

In a bootstrap procedure, the Kármán equation (68) is solved numerically for the initial distribution $R\tau 1(\xi )$ using the *Mathematica*^{TM} function *NDSolve*,

The initial value of *ξ* is chosen to be close to the airfoil leading edge. For the J0012 case, with $U(\xi )$ known analytically, *ξ _{initial}* was chosen to be $10\u22127$. For the NACA 0012 case discussed later, with $U(\xi )$ known approximately,

*ξ*was chosen to be 0.1. The initial distribution, $R\tau 1(\xi )$, is the black curve in Fig. 25 with $R\tau 1TE=9784$.

_{initial}The thickness functions $F0(R\tau 1),\u2009F1(R\tau 1),\u2009F2(R\tau 1)$, and $F3(R\tau 1)$ are computed over the range $0<R\tau 1<R\tau 1TE$ and are used to calculate the first nonzero distribution of the modified Clauser parameter, $\beta c(R\tau 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(\beta c)$ and $b(\beta c)$ are allowed to vary over the airfoil surface according to the correlations (45) and (47).

### B. Step 2

The monotonically increasing friction Reynolds number distribution over the airfoil computed in the first iteration, $R\tau 1(\xi )$ together with $\beta c(R\tau 1)$ are then used to prepare for iteration 2. New thickness functions $F01(R\tau 1),\u2009F11(R\tau 1),\u2009F21(R\tau 1)$, and $F31(R\tau 1)$ are computed over the range $0<R\tau 1<R\tau 1TE$. The Kármán equation (69)

is integrated to generate $R\tau 2(\xi )$. This is the red curve in Fig. 25.

### C. Step 3—Finish

The process is repeated until the *n*th iterate when further changes in $R\tau 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, *F*_{1} and *F*_{2} 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\tau 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\tau iTE(cr)$ at the end of which, $R\tau 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 $\beta c=18$. However, the user needs to be aware that during each iteration, the $\beta ci$ calculated from the current $R\tau i(\xi )$ distribution becomes very large, far larger than 18, as the trailing edge is approached where $u\tau $ 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\tau TE$

*with each iteration*. It turns out that when $\beta ci$ is used to calculate $R\tau i+1$, the calculation reaches the trailing edge at a value of $R\tau i+1TE$ that is considerably less than $R\tau iTE$ with a corresponding $\beta 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=\delta 1/\delta 2$, shown in Figs. 32 and 33. The overall thickness is determined from

## VII. HIGH REYNOLDS NUMBER

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, $\varphi $ Eq. (23), for a range of values of

*β*. Changes in $\varphi $ 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/\delta h$ falling outside of the range $0\u2264y/\delta h\u22641$. The family of shape functions shown in Fig. 36 is accurately approximated by the eighth-order polynomial

_{c}The coefficients $c0(\beta c)$ to $c8(\beta c)$ are approximated by tenth-order polynomial functions of the pressure gradient parameter *β _{c}* and are strictly limited to the range $\u22121.0<\beta c<18.0$ indicated in Fig. 36. The expressions for

*c*

_{0}to

*c*

_{8}are provided in Appendix A. It is important to note that $\varphi (\beta c,y/\delta h)$ in Eq. (71) is universal in the sense that it applies to any pressure gradient wall flow that falls in the range $\u22121.0<\beta c<18.0$. Higher order fits could be used to increase the range of

*β*beyond 18.0.

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

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 $\varphi (ka,m,\beta c,y/\delta 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(\beta c,R\tau ),\u2009F1(\beta c,R\tau ),\u2009F2(\beta c,R\tau )$, and $F3(\beta c,R\tau )$. 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.

R
. _{chord} | $CdvJ0012$ . | $CdvNACA0012UVPintegral$ . | $CdvNACA0012UVPexplicit$ . | $CdvNACA0012SU2$ . | $CdpNACA0012SU2$ . | $CdvLadson80\u2212180\u2009\u2009\u2009\u2009grit$ . |
---|---|---|---|---|---|---|

10^{5} | 0.014 797 0 | 0.014 817 4 | 0.013 499 8 | 0.012 746 0 | 0.004 738 7 | ⋯ |

$5\xd7105$ | 0.010 352 9 | 0.010 397 7 | 0.010 032 5 | 0.010 108 7 | 0.002 664 7 | ⋯ |

10^{6} | 0.009 062 1 | 0.009 147 5 | 0.009 003 4 | 0.008 930 0 | 0.002 252 9 | ⋯ |

$2\xd7106$ | 0.008 040 3 | 0.008 147 7 | 0.008 057 8 | ⋯ | ⋯ | 0.008 53 |

$4\xd7106$ | 0.007 214 1 | 0.007 295 5 | 0.007 244 0 | ⋯ | ⋯ | 0.007 33 |

$5\xd7106$ | 0.006 977 5 | 0.007 050 9 | 0.007 005 5 | 0.007 039 7 | 0.001 608 7 | ⋯ |

$6\xd7106$ | 0.006 787 1 | 0.006 862 6 | 0.006 820 3 | ⋯ | ⋯ | 0.006 82 |

$8.95\xd7106$ | 0.006 397 9 | 0.006 488 3 | 0.006 435 3 | ⋯ | ⋯ | 0.006 51 |

10^{7} | 0.006 296 5 | 0.006 394 3 | 0.006 294 3 | 0.006 431 9 | 0.001 495 7 | ⋯ |

$1.2\xd7107$ | 0.006 137 3 | 0.006 228 17 | 0.006 170 8 | ⋯ | ⋯ | 0.006 53 |

$5\xd7107$ | 0.005 112 8 | 0.005 102 1 | 0.005 082 9 | 0.005 128 0 | 0.001 092 0 | ⋯ |

10^{8} | 0.004 716 9 | 0.004 716 8 | 0.004 6502 | 0.004 696 2 | 0.000 933 7 | ⋯ |

10^{9} | 0.003 512 5 | 0.003 547 7 | 0.003 535 7 | 0.003 400 6 | 0.000 560 9 | ⋯ |

10^{10} | 0.002 755 5 | 0.002 814 7 | 0.002 767 3 | ⋯ | ⋯ | ⋯ |

10^{11} | 0.002 204 9 | 0.002 147 2 | 0.002 217 3 | ⋯ | ⋯ | ⋯ |

10^{12} | 0.001 802 8 | 0.001 764 5 | 0.001 812 6 | ⋯ | ⋯ | ⋯ |

R
. _{chord} | $CdvJ0012$ . | $CdvNACA0012UVPintegral$ . | $CdvNACA0012UVPexplicit$ . | $CdvNACA0012SU2$ . | $CdpNACA0012SU2$ . | $CdvLadson80\u2212180\u2009\u2009\u2009\u2009grit$ . |
---|---|---|---|---|---|---|

10^{5} | 0.014 797 0 | 0.014 817 4 | 0.013 499 8 | 0.012 746 0 | 0.004 738 7 | ⋯ |

$5\xd7105$ | 0.010 352 9 | 0.010 397 7 | 0.010 032 5 | 0.010 108 7 | 0.002 664 7 | ⋯ |

10^{6} | 0.009 062 1 | 0.009 147 5 | 0.009 003 4 | 0.008 930 0 | 0.002 252 9 | ⋯ |

$2\xd7106$ | 0.008 040 3 | 0.008 147 7 | 0.008 057 8 | ⋯ | ⋯ | 0.008 53 |

$4\xd7106$ | 0.007 214 1 | 0.007 295 5 | 0.007 244 0 | ⋯ | ⋯ | 0.007 33 |

$5\xd7106$ | 0.006 977 5 | 0.007 050 9 | 0.007 005 5 | 0.007 039 7 | 0.001 608 7 | ⋯ |

$6\xd7106$ | 0.006 787 1 | 0.006 862 6 | 0.006 820 3 | ⋯ | ⋯ | 0.006 82 |

$8.95\xd7106$ | 0.006 397 9 | 0.006 488 3 | 0.006 435 3 | ⋯ | ⋯ | 0.006 51 |

10^{7} | 0.006 296 5 | 0.006 394 3 | 0.006 294 3 | 0.006 431 9 | 0.001 495 7 | ⋯ |

$1.2\xd7107$ | 0.006 137 3 | 0.006 228 17 | 0.006 170 8 | ⋯ | ⋯ | 0.006 53 |

$5\xd7107$ | 0.005 112 8 | 0.005 102 1 | 0.005 082 9 | 0.005 128 0 | 0.001 092 0 | ⋯ |

10^{8} | 0.004 716 9 | 0.004 716 8 | 0.004 6502 | 0.004 696 2 | 0.000 933 7 | ⋯ |

10^{9} | 0.003 512 5 | 0.003 547 7 | 0.003 535 7 | 0.003 400 6 | 0.000 560 9 | ⋯ |

10^{10} | 0.002 755 5 | 0.002 814 7 | 0.002 767 3 | ⋯ | ⋯ | ⋯ |

10^{11} | 0.002 204 9 | 0.002 147 2 | 0.002 217 3 | ⋯ | ⋯ | ⋯ |

10^{12} | 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$.^{17} Figure 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\tau $, as a function of $R\tau $. Above $R\tau =105$ the error is less than $0.5%$.

## VIII. VISCOUS DRAG COEFFICIENT OF THE JOUKOWSKY AND NACA 0012 AIRFOILS

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

with constants, $a1=0.177\u2009349\u2009856,\u2009a2=\u22120.075\u2009600\u2009000,\u2009a3=\u22120.212\u2009843\u2009959\u20091,\u2009a4=0.173\u2009640\u2009303\u20090$, and $a5=\u22120.062\u2009546\u2009200\u20092$, 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(\xi \u0303)$ determined using conformal mapping, compared with the free stream $U(\xi \u0303)$ computed from the SU2 pressure distribution at $Rchord=107$. The pressure drag from SU2 at this Reynolds number was $CdpSU2=0.001\u2009496$.

The final *C _{dv}* vs

*R*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\xaf,a\xaf,m\xaf)$, 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

_{chord}*C*at $R\tau =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

_{f}*σ*'s given in Table I. The position of the magenta error bar at $Rchord=7\xd7106$ is approximately at $R\tau =5000$. It should be noted that, while the error associated with changing the optimal parameters this way actually depends on $R\tau $, 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.

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 Ladson^{40} who used carborundum strips placed at $x\u0303/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 McCroskey^{41}

Note that an interpolation of the pressure drag coefficient, $CdpSU2$ is subtracted from the *C _{d}* data in Ladson

^{40}and McCroskey

^{41}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 SU2^{33} using the Spalart–Allmaras model.^{37} No transition model was used in the computations. Between $Rchord=106$ and 10^{9}, 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 *R _{chord}* cases. The finest grid available from the Langley Turbulence Modeling Resource webpage

^{42}gave a wall spacing of $10\u22127$ m corresponding to a minimum $y+\u22731$ over significant portions of the airfoil instead of the preferred $y+\u22641$. 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.

### A. Drag force comparison between the UVP and SU2 NACA0012 results at $Rchord=107$

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.

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.

## IX. CONCLUDING REMARKS

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

*C*with

_{dv}*b*and

*n*fixed at the zero pressure gradient values is remarkably close to the final value, although the distributions of $R\tau $ and

*C*generated by the modified Clauser parameter are essential to properly capture the trailing edge flow.

_{f}For all calculations, the parameters that characterize the wall flow, $(k\xaf,a\xaf,m\xaf)$ 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 *e ^{N}*.

^{43,44}For a boundary layer developing in a very low free stream turbulence environment, for example, Schubauer and Klebanoff,

^{30}

*N*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\u22641$. 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.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX A: SHAPE FUNCTION COEFFICIENTS

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

### APPENDIX B: UVP INTEGRAL METHOD FLOW CHART

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