This work explores simple relations that follow from the momentum-integral equation absent a pressure gradient. The resulting expressions enable us to relate the boundary-layer characteristics of a velocity profile, *u*(*y*), to an assumed flow function and its wall derivative relative to the wall-normal coordinate, *y*. Consequently, disturbance, displacement, and momentum thicknesses, as well as skin friction and drag coefficients, which are typically evaluated and tabulated in classical monographs, can be readily determined for a given profile, *F*(*ξ*) = *u*/*U*. Here, *ξ* = *y*/*δ* denotes the boundary-layer coordinate. These expressions are then employed to provide a rational explanation for the 1921 Pohlhausen polynomial paradox, namely, the reason why a quartic representation of the velocity leads to less accurate predictions of the disturbance, displacement, and momentum thicknesses than using cubic or quadratic polynomials. Not only do we identify the factors underlying this behavior but also we proceed to outline a procedure to overcome its manifestation at any order. This enables us to derive optimal piecewise approximations that do not suffer from the particular limitations affecting Pohlhausen’s *F* = 2*ξ* − 2*ξ*^{3} + *ξ*^{4}. For example, our alternative profile, *F* = (5*ξ* − 3*ξ*^{3} + *ξ*^{4})/3, leads to an order-of-magnitude improvement in precision when incorporated into the Kármán–Pohlhausen approach in both viscous and thermal analyses. Then, noting the significance of the Blasius constant, $s\xaf\u22481.630\u2009398$, this approach is extended to construct a set of uniformly valid solutions, including $F=1\u2212exp[\u2212s\xaf\xi (1+12s\xaf\xi +\xi 2)]$, which continues to hold beyond the boundary-layer edge as *y* → ∞. Given its substantially reduced error, the latter is shown, through comparisons to other models, to be practically equivalent to the Blasius solution.

## I. INTRODUCTION

Following Prandtl’s pioneering treatment of boundary layers^{1} and the shape-preserving similarity solution reported by Blasius,^{2} the Kármán–Pohlhausen (KP) momentum-integral approach, introduced through two sequential papers by Kármán^{3} and Pohlhausen,^{4} is widely regarded as the staple in boundary-layer analysis. This is mainly due to its effectiveness at providing a wealth of useful details for boundary layers in both low and high Reynolds number flows. Despite its simplicity, the KP approach continues to play a pivotal role in several monographs that devote themselves to boundary-layer theory. These include those by Rosenhead,^{5} Schlichting,^{6} Cebeci and Cousteix,^{7} Oleinik and Samokhin,^{8} White,^{9} Schetz and Bowersox,^{10} Pritchard and Mitchell,^{11} Schlichting and Gersten,^{12} and others. It must be noted, however, that the KP approach is restricted to the analysis of the so-called inner region, with the far field referring to the outer inviscid freestream.

Given its ease of implementation and broad-brush bookkeeping outcomes in a variety of viscous-flow problems, the KP approach has been applied to a plethora of thermofluid flow configurations. Examples abound and one may cite studies that involve extensions to Ostwald–de Waele power-law fluids,^{13} rarefied gases over porous plates,^{14} rotating three-dimensional boundary layers on cones,^{15} thermal boundary layers on short rotating blades,^{16} forced convection from surfaces immersed in non-Newtonian fluids using planar,^{17} axisymmetric,^{18} and irregular geometry,^{19} forced convection in highly pseudoplastic fluids,^{20} natural convection over bodies immersed in fluid-saturated porous media,^{21} non-Newtonian motions within slider bearings,^{22} natural convection over bodies embedded in porous media,^{23} non-Newtonian boundary layers for low-Reynolds number mud flows,^{24} radial film flows,^{25} asymmetrically mixed convective motions around horizontal cylinders,^{26} heat transfer from infinite isothermal cylinders with both Newtonian and non-Newtonian fluids,^{27,28} momentum and thermal analyses of power-law fluids over flat plates,^{29} and forced convection in power-law fluid-saturated porous media.^{30}

The KP approach has also been employed at the basis of predicting the stability of laminar boundary layers at separation,^{31} although solutions of the differential equations of motion, often manifested in the similarity-preserving Blasius^{32–39} and Falkner–Skan type expressions,^{40–42} have gradually superseded integral approximations. The resulting methodological shift may be attributed, in part, to the increased precision often associated with the differential approach and, in part, to the substantial advancements made in computational tools.

In the present investigation, however, the objective is not to extend the KP approach to a new flow configuration. Rather, it stands in reducing the KP formulation in the absence of a pressure gradient to simple relations that display direct connections between their predictions and their true values for flat-plate motion. Our next objective is to revisit the assumptions made by Pohlhausen,^{4} which have been used at the basis of constructing closed-form analytic expressions for a series of conjectured inner velocity profiles; the latter may be separated into two principal categories: those seeking piecewise approximations, valid only within the boundary layer,^{4,43} and those pursuing monotonically increasing continuous profiles that merge smoothly into the far field as *y* → ∞.^{44–51} We begin by considering those that extend over 0 ≤ *y* ≤ *δ*, where *y* = *δ* represents the edge of the viscous layer and *y* denotes the wall-normal coordinate. As usual, we assign *U* to the freestream velocity and provide a labeled sketch of the streamline-bounded control volume and coordinate system in Fig. 1.

Since the majority of the continuous flow formulations aim at producing uniformly valid analytical solutions to the Blasius equation, we find it necessary to evaluate the viability of the corresponding boundary conditions, which are embedded in a KP formulation, *vis-à-vis* those predicted by a differential technique, to retrieve comparable velocity approximations. This two-pronged effort is carried out while taking into account the subtle connection between the Blasius similarity variable, *η*, and the normalized coordinate, *ξ*, which arises in the integral formulation. By closely examining parallel predictions from the integral and differential approaches for the same flow configuration, we are able to identify previously used boundary conditions that, when replaced by more suitable constraints, lead to substantially more accurate flow profiles. More specifically, we show that Pohlhausen’s quartic polynomial, *F* = *u*/*U* = 2*ξ* − 2*ξ*^{3} + *ξ*^{4}, which has been routinely used in viscous and thermal analyses, with and without pressure gradients, can be definitively amended by allowing it to satisfy a more realistic boundary condition at the edge of the boundary layer. By doing so, we arrive at an alternative representation, *F* = (5*ξ* − 3*ξ*^{3} + *ξ*^{4})/3, which will be shown to outperform other piecewise approximations in its prediction of both viscous and thermal boundary-layer characteristics.

Recognizing the utility of this approach in deriving spatially restricted piecewise approximations and since no flat-plate study is complete with no referral to the Blasius equation, we circle back and apply the same procedure to obtain a continuous approximation to the Blasius problem over its semi-infinite domain. The latter continues to serve as an important benchmark and testing ground for new viscous and heat transfer approximations as well as stability and transition analyses entailing both traditional and non-Newtonian fluids, with some including wall roughness, inclination, and suction or injection.^{52–62}

By taking advantage of the Blasius constant at the wall, $F\u2032(0)=s\xaf\u22481.630\u2009398$, and the crucial role that it plays in prescribing the trajectory of the solution under construction, we identify five essential boundary conditions that are inspired by Pohlhausen’s classic requirements; the latter are sequentially applied in a manner to produce a compact, decaying exponential function that entails a maximum *L*_{2} error of 6.4 × 10^{−4}. Such a small discrepancy places it within the same margin of error affecting the Blasius equation itself. We close by extending this approach to higher orders and by comparing the one-term expression, $F=1\u2212exp[\u2212s\xaf\xi (1+12s\xaf\xi +\xi 2)]$, to other available profiles in the semi-infinite range of 0 ≤ *ξ* < ∞. Our verification tools and error analyses are systematically used to demonstrate that the decaying exponential formulation matches the Blasius solution virtually identically, not only from the wall to the edge of the boundary layer but also all the way to infinity.

From an organizational standpoint, this paper is divided into three main sections, besides the Introduction and Conclusions. In Sec. II, several reduced momentum-integral relations are presented to help connect the fundamental boundary-layer properties, including the wall-normal velocity and vorticity, directly to an assumed velocity profile. This is followed in Sec. III with a procedure to rapidly evaluate boundary-layer properties using piecewise velocity approximations. Along the way, the paradoxical behavior of Pohlhausen’s high-order polynomial approximations is brought to light and discussed. The specific factors leading to unexpected deviations from the correct Blasius predictions with increasing polynomial orders are also identified. By introducing an alternative set of boundary conditions, the ability to construct high-order polynomials with increasing degrees of accuracy is demonstrated, starting with a quartic polynomial. The latter is shown to reduce the overall *L*_{2} error associated with the KP method by one order of magnitude. Next, the improved effectiveness of the KP approach, when used in concert with more realistic velocity and temperature profiles, is confirmed in the context of thermal analysis. More specifically, in the case of flat-plate heat transfer, it is shown that the KP approach is capable of predicting the local Nusselt number within a half percent relative error. Finally, the same tools developed in Secs. II and III are extended in Sec. IV to the systematic development of continuous decaying exponential profiles that mirror the characteristic behavior of the Blasius solution almost flawlessly. The accuracy of these simple approximations relative to other comparable models in the literature is also examined and discussed.

## II. REDUCED MOMENTUM-INTEGRAL RELATIONS

In this section, the geometry and KP formulation are presented in the context of flat-plate analysis in conjunction with the use of a conjectured viscous-flow approximation for *u*(*y*). Direct relations that link *u*(*y*), or more generally, *F*(*ξ*), to traditionally computed properties, such as the disturbance, displacement, and momentum thicknesses, are also established.

### A. Flat-plate configuration

Our schematic of the physical setting in Fig. 1 sets the sharp edge of the plate at (*x*, *y*) = (0, 0) and shows the fluid shearing at the wall as it sweeps along the plate, where a tangential stress, *τ*_{$w$}, is exerted. The velocity distribution, *u*(*y*), at any downstream station, *x*, exhibits a smooth drop-off before vanishing at the wall. The normal “upwelling” velocity is taken to be $v$, and the disturbance or shear-layer thickness is labeled *δ*. Streamlines situated outside this layer are deflected vertically by a distance *δ*^{*}, namely, the displacement thickness. The upper streamline slings upward from *y* = *h* at *x* = 0 to *y* = *δ* = *h* + *δ*^{*} at *x* = *L* in a control volume of length *L*. In this work, we assume that the freestream just upstream of the sharp edge is given by *u* = *U* = const. As usual, we normalize the velocity and transverse coordinate using *F* ≡ *u*/*U* and *ξ* ≡ *y*/*δ*.

### B. Reduced momentum-integral relations

By focusing on the incompressible planar motion over a flat plate at zero incidence, one can use continuity to determine $h=\u222b0\delta (u/U)dy=\delta \u222b01F(\xi )d\xi $. Then, using *ρ* for density, a momentum balance enables us to write, as detailed by Schlichting^{6} or White,^{9}

where *C*_{f} stands for the skin friction coefficient. The incompressible displacement and momentum thicknesses can be expressed as

We recall that the semi-infinite to finite domain approximations are permitted because of the negligible contributions of the integrands for *y* > *δ*. In fact, these approximations become strict identities for piecewise velocity profiles that remain equal to unity for *y* ≥ *δ*. Moreover, these expressions hold true for any incompressible fluid, be it laminar or turbulent, with constant or variable pressure. We also find it useful to define the non-dimensional displacement and momentum thicknesses, *α* and *β*, along with the shape factor, *H*, using

and

where $\alpha \u0303$ and $\beta \u0303$ represent the approximated integral values of *α* and *β* over the finite [0, 1] interval. Note the strict equalities of $\alpha \u0303=\alpha $ and $\beta \u0303=\beta $ for piecewise analytic profiles exhibiting the property *F*(*ξ* ≥ 1) = 1. For this reason, the distinction of using $\alpha \u0303$ or $\beta \u0303$ is only necessary in the calculation of integral properties associated with uniformly valid, continuously analytic profiles for which the default integrals are evaluated over a semi-infinite domain. For profiles that continue to increase in the far field when *ξ* > 1, restoring the semi-infinite bounds on *α* and *β* in Eqs. (2) and (3) is necessary to avoid underestimating their values. In the case of the Blasius solution, for example, the semi-infinite bounds will result in 0.284% and 0.733% higher values for *α* and *β*, respectively. We also recall that *α*, *β*, and *H* are pure constants that depend strictly on the velocity distribution *F*. Since the focus of this study is on the steady form of Eq. (1), we continue with

which, for a spatially invariant freestream, becomes

### C. Direct connections to a conjectured velocity profile

To make further headway, the following procedure may be used. First, a piecewise analytic velocity may be posited in the form

Once *F* is specified, the constants *α* and *β* may be readily determined from their simple integrals in Eq. (4). Then, using the stress–strain relation at the wall, *τ*_{$w$} can be written as a function of the disturbance thickness and the gradient of the assumed velocity at *y* = *ξ* = 0. On the one hand, we recover

Here, *μ* is the dynamic viscosity and primes indicate differentiation with respect to *ξ*. On the other hand, *τ*_{w} may be related to the axial gradient of the disturbance thickness based on Eq. (6),

These two expressions may be equated to obtain

with *ν* being the kinematic viscosity. For a boundary layer that starts to develop at *x* = 0, one can use *δ*(0) = 0 and integrate Eq. (10) from the plate’s leading edge to any station *x*. After some rearranging, one recovers the traditional dependence on the local Reynolds number, namely,

where *Re*_{x} = *Ux*/*ν*. Note the clear connection between the traditional constant *a* and both the non-dimensional momentum thickness, *β*, and the initial slope, *s* = *F*′(0). For the entire class of polynomial approximations of the form *F* = *a*_{0} + *a*_{1}*ξ* + *a*_{2}*ξ*^{2} + ⋯, *s* is simply *a*_{1}, which can be found by inspection. Another substitution into the friction coefficient in Eq. (6) leads to

The remaining boundary-layer properties follow suit. Recalling that the displacement and momentum thicknesses are deducible from their reversed definitions, *δ*^{*} = *δα* and *θ* = *δβ*, we get

This outcome confirms the equivalence of *θ*/*x* and the skin friction coefficient in Eq. (12). Interestingly, we find the products and ratios of the traditional constants *b* and *a*, which are often tabulated in various monographs,^{9–12} to be proportional to *s* and *β*,

It follows that

Practically, for a given profile *F*, it is only necessary to calculate *α*, *β*, and *s* to deduce $a=2s/\beta $, *b* = *d* = *aβ*, and *c* = *aα*. As for the global drag coefficient, it may be evaluated using^{3}

In view of Eq. (12), it is clear that *C*_{D} = 2*C*_{f}.

At this point, we can revert back to the continuity equation for the purpose of retrieving valuable estimates for $v$ and its maximum value at the boundary-layer edge. We start by switching *u*(*y*) to *F*(*ξ*) in the continuity integral,

This enables us to identify the presence of *δ*(*x*), which may be readily expanded and rearranged within the integrand using

Moreover, by virtue of Eq. (11), the axial gradient of *δ*(*x*) can be eliminated in favor of the local Reynolds number. The resulting expression can be further normalized by *U* and integrated by parts. One gets

where the normal velocity’s characteristic function *G*(*ξ*) can be determined directly from

Naturally, *G*(*ξ*) will depend on the profile at hand. For piecewise velocity representations with *F*(*ξ* ≥ 1) = 1, the peak normal velocity may be shown to occur at *y* = *δ*, where

To determine *G*_{δ} explicitly, it is helpful to recreate the definition of $\alpha \u0303$ in Eq. (20). By simple manipulation, we arrive at

This expression is somewhat illuminating. It shows that for all assumed velocity approximations that observe the *F*(*ξ* ≥ 1) = 1 condition, the maximum normal velocity coefficient reduces to $G\delta =a\alpha \u0303/2=a\alpha /2=c/2$ by way of Eq. (13). As such and in reference to standard skin friction and momentum thickness tables,^{9–12} the peak normal velocity coefficients, $Rex(v/U)max$, are simply half of the tabulated values for $\delta *Rex/x$.

However, for those approximations that observe Prandtl’s *F*(1) = 0.99 cutoff requirement on *u*/*U*, we get $G\delta =a\alpha \u0303/2\u22120.005a$. Moreover, for uniformly valid velocity profiles that continue to increase past *ξ* = 1, the maximum normal velocity occurs at infinity rather than the edge of the disturbance thickness. Besides the normal velocity at *ξ* = 1, a slightly larger value may be realized at

Here, too, the maximum *G*_{∞} coefficient may be readily deduced from Eq. (20). Making use of *F*(∞) = 1, we can add and subtract the variable *ξ* by putting

This useful outcome is yet contingent upon $\xi F(\xi )\u22121\u21920$ as *ξ* → ∞. Such will be the case when

The last condition in Eq. (25) is safely secured because, in boundary-layer theory, the gradient of the velocity in the far field is known to approach zero exponentially, i.e., faster than any polynomial order. Furthermore, the use of ($v$/*U*)_{max} = *aα*/2 is consistent with corresponding estimates reported in the literature.^{9–12} Another benefit of these relations is that they enable us to link the key characteristic properties explicitly to *a*, *α*, and *F*.

At this juncture, having determined both *u* and $v$, we can proceed by evaluating the boundary-layer vorticity, *ω*. Being a cross-derivative function, the outcome is straightforward. We get

or

where *ω*^{*} denotes the non-dimensional vorticity. The last expression may be further simplified using Eq. (20) into

Compared to the normal velocity and its peak value in Eqs. (19) and (21), which are inversely proportional to the square root of the Reynolds number, the *G*′-dependent term in Eq. (26) is appreciably smaller, being proportional to $Rex\u22121$. As such, *ω*^{*} ≈ −*F*′(*ξ*) is a suitable one-term approximation for the vorticity, whose peak value at the wall, *ω*^{*}(0) = −*s*, becomes the negative of the so-called “connection parameter” in the corresponding Blasius problem. This equality is consistent with the foundational role that shear-layer vorticity plays in prescribing the behavior of fluid motion.

Before leaving this section, it may be helpful to reiterate that the procedure just described is fairly straightforward and facilitates the evaluation of most salient boundary-layer features for a given velocity distribution. Once the defining integrals for *α* and *β* are determined and *s* is in hand, the remaining properties can be calculated with minimal effort using Eqs. (11)–(22), (26), and (27). In fact, the ability to connect various boundary-layer properties directly to an assumed velocity profile will prove essential to the remaining segments of this work.

## III. IMPROVED EFFECTIVENESS OF THE MOMENTUM-INTEGRAL APPROACH

For decades, integral methods have been dominated by the use of intuitive, guessed velocity profiles. As indicated earlier, the Kármán–Pohlhausen approach is named as such after two concurrent articles by Kármán,^{3} who provided the integral formulation, and Pohlhausen,^{4} who presented a systematic procedure for constructing and integrating suitable velocity approximations. A side benefit of the KP approach is that once a reasonable form of *u* is conjectured, accurate estimates of the boundary-layer properties can be expected because integration tends to wash out the positive and negative deviations of the assumed profile from the exact velocity distribution. However, this procedure is not without limitations.

### A. Direct evaluation of properties using piecewise velocity approximations

To illustrate the utility of the foregoing relations, the boundary-layer properties will be evaluated using several closed-form approximations that have been widely used in the literature. The accuracy of these profiles will then be gauged relative to the classic Blasius solution for flow over a flat plate. We begin by examining five piecewise models that are defined for 0 ≤ *y* ≤ *δ* and that are set equal to unity over the semi-infinite domain, *δ* < *y* < ∞. These correspond to Pohlhausen’s quadratic, cubic, and quartic profiles,^{4}

Schlichting’s sinusoidal profile,^{43}

and a rationally optimized quartic profile,

Using their normalized forms, these profiles are collected in Table I and compared to the classic solution due to Blasius.^{2} Note that relative percentage errors are provided below each estimate along with the overall $L2=[\u222b01(F\u2212F\xaf)2d\xi ]1/2$ error across the viscous layer, where a “barred” quantity refers to its Blasius value. This is performed to help quantify the overall deviation in each approximation from the traditionally reported Blasius values.^{11,12} Properties that are not shown can be easily deduced because, for laminar flow, $(\theta /x)Rex=CfRex=b$, which is already tabulated, and $CDReL$ is twice this constant, as per Eq. (16). Conversely, the maximum wall-normal velocity, $Rex(v/U)max$, is simply half of the entries for $\delta *Rex/x$, as per Eq. (24).

$F(\xi )=uU$ . | $\alpha =\delta *\delta $ . | $\beta =\theta \delta $ . | $H=\delta *\theta $ . | $\delta xRex$ . | $CfRex$ . | $\delta *xRex$ . | L_{2} error
. |
---|---|---|---|---|---|---|---|

2ξ − ξ^{2} | 0.333 | 0.133 | 2.500 | 5.477 | 0.730 | 1.826 | 0.020 |

3.1% | 0.25% | 3.5% | 9.5% | 10% | 6.1% | ||

$32\xi \u221212\xi 3$ | 0.375 | 0.139 | 2.692 | 4.641 | 0.646 | 1.740 | 0.034 |

9.0% | 4.7% | 4.0% | 7.2% | 2.6% | 1.1% | ||

2ξ − 2ξ^{3} + ξ^{4} | 0.300 | 0.118 | 2.554 | 5.836 | 0.685 | 1.751 | 0.054 |

13% | 12% | 1.4% | 17% | 3.2% | 1.8% | ||

$sin12\pi \xi $ | 0.363 | 0.137 | 2.660 | 4.795 | 0.655 | 1.743 | 0.021 |

5.6% | 2.7% | 2.7% | 4.1% | 1.3% | 1.3% | ||

$53\xi \u2212\xi 3+13\xi 4$ | 0.350 | 0.134 | 2.618 | 4.993 | 0.668 | 1.748 | 0.008 |

1.7% | 0.52% | 1.1% | 0.13% | 0.53% | 1.6% | ||

Blasius^{2} | 0.344 | 0.133 | 2.59 | 5 | 0.664 | 1.72 |

$F(\xi )=uU$ . | $\alpha =\delta *\delta $ . | $\beta =\theta \delta $ . | $H=\delta *\theta $ . | $\delta xRex$ . | $CfRex$ . | $\delta *xRex$ . | L_{2} error
. |
---|---|---|---|---|---|---|---|

2ξ − ξ^{2} | 0.333 | 0.133 | 2.500 | 5.477 | 0.730 | 1.826 | 0.020 |

3.1% | 0.25% | 3.5% | 9.5% | 10% | 6.1% | ||

$32\xi \u221212\xi 3$ | 0.375 | 0.139 | 2.692 | 4.641 | 0.646 | 1.740 | 0.034 |

9.0% | 4.7% | 4.0% | 7.2% | 2.6% | 1.1% | ||

2ξ − 2ξ^{3} + ξ^{4} | 0.300 | 0.118 | 2.554 | 5.836 | 0.685 | 1.751 | 0.054 |

13% | 12% | 1.4% | 17% | 3.2% | 1.8% | ||

$sin12\pi \xi $ | 0.363 | 0.137 | 2.660 | 4.795 | 0.655 | 1.743 | 0.021 |

5.6% | 2.7% | 2.7% | 4.1% | 1.3% | 1.3% | ||

$53\xi \u2212\xi 3+13\xi 4$ | 0.350 | 0.134 | 2.618 | 4.993 | 0.668 | 1.748 | 0.008 |

1.7% | 0.52% | 1.1% | 0.13% | 0.53% | 1.6% | ||

Blasius^{2} | 0.344 | 0.133 | 2.59 | 5 | 0.664 | 1.72 |

Apart from the optimized quartic profile, whose error does not exceed 1.7% in the 0 ≤ *ξ* ≤ 1 interval, it may be seen that the maximum error in each case varies between 6% and 17%, which is typical of integral theories. Some profiles, however, appear to be distinctly more precise than others, namely, in mirroring the behavior of the Blasius solution. What is most perplexing, perhaps, is what some academicians have come to identify, rather informally, as the Pohlhausen paradox: one may notice that as we move from Pohlhausen’s second-order to his fourth-order polynomial—which is capable of securing additional boundary conditions—the agreement with the Blasius solution does not improve, as one would expect, but rather deteriorates. Being part of a widely taught subject, this inconsistency has often puzzled fluid mechanics learners and instructors alike. From Table I, it can be seen that the errors in estimating the non-dimensional displacement and momentum thicknesses {*α*, *β*} leap by one order of magnitude, i.e., from {3.1%, 0.25%} to {13%, 12%}, while the overall *L*_{2} error increases from 0.02 to 0.05, when Pohlhausen’s quadratic polynomial is exchanged with its quartic counterpart. The error in predicting the boundary-layer thickness also increases from 9.5% to a whopping 17%.

### B. Paradoxical behavior of Pohlhausen’s polynomials

Hence, one may wonder what is controlling the accuracy of these profiles? For nearly a century, the underlying assumption has been that the fidelity of the conjectured profiles could be improved by securing additional boundary conditions that are representative of the Blasius model. In this vein, Pohlhausen^{4} proposed five cardinal requirements that, apart from the wall adherence and smooth merging with the freestream conditions, i.e., the three constraints satisfied by the quadratic polynomial in Eq. (28), piled on two further requirements on the shear stress: one at the wall and the other one at the edge of the viscous domain. To summarize, the first three foundational constraints consist of

and

As for the fourth condition, it serves to maintain the correct momentum balance at the wall between the non-vanishing components of the boundary-layer equation via

Here, the pressure gradient within the boundary layer is expressed in terms of the freestream velocity using Euler’s equation. One may verify that not only do Pohlhausen’s cubic and quartic polynomials satisfy this momentum balance requirement but also do Schlichting’s sinusoidal and the optimized quartic profiles.

Finally, Pohlhausen’s fifth condition, i.e., the one that we scrutinize here, assumes zero shear at the edge, which translates into a vanishing curvature at *y* = *δ*,

Of the five profiles described so far, only Pohlhausen’s quartic polynomial satisfies this condition identically. Incidentally, it also proves to be the most imprecise, as it leads to the largest overall *L*_{2} error and a 17% discrepancy in predicting the boundary-layer thickness relative to the traditionally reported Blasius values.^{9–12}

As it turns out, although the fifth condition holds true as *y* → ∞, it proves to be unsuitable at the edge of the boundary layer where the velocity gradient continues to change. In fact, a sufficiently resolved numerical solution of the Blasius equation with respect to the boundary-layer coordinate *ξ* returns a curvature of −0.7085, which is not yet zero, but rather of order unity (Fig. 2).

Subjecting the quartic polynomial to the fifth constraint, which happens to be off by one order of magnitude, partly explains its tendency to overpredict several boundary-layer properties relative to its lower-order models (cf. Table I). It also leads to a higher initial slope at *y* = *ξ* = 0, as one may readily infer from Fig. 2(a); therein, the piecewise analytic profiles are superimposed on the Blasius line. In contrast, Schlichting’s sinusoidal and the optimized quartic profiles display slopes that are closer to the exact Blasius value of 1.630 at *y* = 0. This may be attributed to their wall gradients of *F*′(0) = *π*/2 ≈ 1.571 and 5/3 ≈ 1.667, which slightly undershoot and overshoot the true slope, respectively. Other deviations from the Blasius solution can be detected in the two portions of the graph, especially in the magnified inset of Fig. 2(b). Interestingly, the assumption that *u* = *U* (therefore, *F* = 1 instead of 0.99) at *y* = *δ* (far right) can be seen to affect all of the piecewise approximations in Fig. 2(a), except for the Blasius curve. The latter returns a value of 0.99 at *y* = *δ*, which is consistent with the 99% disturbance basis. This discrepancy and how to overcome it will be explored in Sec. IV C.

For further clarity, the properties and conditions observed by the five piecewise profiles and their derivatives at the endpoints of the viscous domain are summarized in Table II. Although not imposed, the negative curvature of the quadratic Pohlhausen profile at *ξ* = 1, which matches the curvature of the optimized quartic polynomial of −2, appears to be more in-line with the computed Blasius edge curvature than its quartic companion. This characteristic feature helps explain, in part, its improved accuracy relative to its quartic form, despite its inability to secure the momentum balance requirement at the wall, i.e., Pohlhausen’s fourth condition [Eq. (34)]. In fact, it is only after recognizing these various limitations that one can manage to derive a rationally optimized quartic polynomial that satisfies Pohlhausen’s four essential boundary conditions while judiciously avoiding the fifth requirement. To make further headway, we posit that the main defect in the classic quartic solution is fundamentally triggered by a prematurely imposed zero shear-stress condition at *y* = *δ*.

F(ξ) = u/U
. | F(0)
. | F′(0)
. | $F\u2032\u2032$(0) . | F(1)
. | F′(1)
. | $F\u2032\u2032$(1) . |
---|---|---|---|---|---|---|

2ξ − ξ^{2} | 0 | 2.000 | −2 | 1.000 | 0 | −2.000 |

$32\xi \u221212\xi 3$ | 0 | 1.500 | 0 | 1.000 | 0 | −3.000 |

2ξ − 2ξ^{3} + ξ^{4} | 0 | 2.000 | 0 | 1.000 | 0 | 0.000 |

$sin(12\pi \xi )$ | 0 | 1.571 | 0 | 1.000 | 0 | −2.467 |

$53\xi \u2212\xi 3+13\xi 4$ | 0 | 1.667 | 0 | 1.000 | 0 | −2.000 |

Blasius^{2} | 0 | 1.630 | 0 | 0.990 | 0.0904 | −0.709 |

F(ξ) = u/U
. | F(0)
. | F′(0)
. | $F\u2032\u2032$(0) . | F(1)
. | F′(1)
. | $F\u2032\u2032$(1) . |
---|---|---|---|---|---|---|

2ξ − ξ^{2} | 0 | 2.000 | −2 | 1.000 | 0 | −2.000 |

$32\xi \u221212\xi 3$ | 0 | 1.500 | 0 | 1.000 | 0 | −3.000 |

2ξ − 2ξ^{3} + ξ^{4} | 0 | 2.000 | 0 | 1.000 | 0 | 0.000 |

$sin(12\pi \xi )$ | 0 | 1.571 | 0 | 1.000 | 0 | −2.467 |

$53\xi \u2212\xi 3+13\xi 4$ | 0 | 1.667 | 0 | 1.000 | 0 | −2.000 |

Blasius^{2} | 0 | 1.630 | 0 | 0.990 | 0.0904 | −0.709 |

### C. Procedure to develop increasingly accurate polynomials

The long-standing paradox associated with Pohlhausen’s polynomials, which are incorporated into the KP approach, is that their accuracy tends to deteriorate at increasing orders. This behavior is manifestly established in Table I, where their consistently increasing *L*_{2} errors with successive increases in their polynomial orders are documented. This is also evidenced by their increasing spatial deviations from the Blasius solution in Fig. 2, including its magnified inset. Although counter-intuitive and thought-provoking paradoxes are somewhat popular in fluid analysis, a formal procedure that overcomes this perplexing behavior is desirable.

To begin, we opt for a non-dimensional velocity profile that can be approximated by an *N*th-order polynomial of the form

Based on Eq. (4), the non-dimensional displacement and momentum thicknesses can be written as

Polynomials that conform to this series expansion and that satisfy some (or all) of the conditions given by Eqs. (31)–(34) include, but are not limited to,

On the one hand, when Pohlhausen’s physical requirements in Eqs. (31)–(34) are applied on *F*, they translate into

On the other hand, a robust numerical solution of the Blasius equation predicts, as shown in Table II,

and

It can thus be seen that the $F\u2032\u2032$(1) = 0 constraint, which is secured by Pohlhausen’s quartic solution, entails a non-negligible error. To overcome this limitation, we hypothesize that the $F\u2032\u2032$(1) = 0 requirement must be released. The degree of freedom thereby gained can be leveraged while taking into account the direct relations identified in Sec. II C, i.e., to produce the closest values of {*α*, *β*, *s*} to the problem’s exact solution.

To illustrate this procedure, we choose a quartic polynomial that satisfies all conditions in Eq. (38) except $F\u2032\u2032$(1) = 0. Despite the penalty that the *F*′(1) = 0 assumption engenders, it is retained in observance of Pohlhausen’s tradition of suppressing the gradient of the velocity at the edge of the domain. Such a profile consists of

Then, using the reduced relations specified in Eq. (4), the non-dimensional disturbance and momentum thicknesses can be expressed explicitly as functions of *s*. We get

At the time the KP approach was being developed,^{63} the classical Blasius solution predicted $\alpha \xaf\u22480.344$, $\beta \xaf\u22480.133$, and $s\xaf\u22481.660$.^{9–12} To minimize the overall error relative to the traditionally reported Blasius values, it is useful to define the following objective function:

We have chosen to constrain the total relative residual of *α*, *β*, and *s* because these parameters have been shown in Sec. II C to be the most influential in prescribing the remaining boundary-layer properties. At this juncture, we may substitute Eq. (41) into Eq. (42) and express $R$ solely in terms of *s*. This enables us to minimize $R(s)$ analytically, namely, by setting

where

In this manner, by differentiating the objective function with respect to *s*, one is rationally guided to select *s*^{*} ≈ 1.6794 ≈ 5/3, which is still 2.2% higher than the Blasius slope of $s\xaf\u22481.6304$ and, unsurprisingly, lower than the classic Pohlhausen slope of 2 by precisely 20%. Evidently, the use of *s*^{*} leads to an excellent overall approximation, albeit slightly different from the Blasius constant. Finally, by inserting this trajectory back into Eq. (40) and simplifying, we arrive at

It is a simple matter to verify that this rationally optimized profile satisfies Pohlhausen’s four fundamental boundary conditions while dodging the unrealistic $F\u2032\u2032$(1) = 0 constraint. One may also confirm that this solution leads to reliable estimates for the boundary-layer thickness and skin friction, namely, $(\delta /x)Rex\u22484.993$ and $CfRex\u22480.668$. These estimates differ from the classical Blasius values of 5.0 and 0.664 by 0.13% and 0.53%, respectively. Compared to other profiles featured in Table I, its *L*_{2} error of 0.008, evaluated over its finite interval, is another favorable metric: it reflects an order of magnitude reduction relative to its predecessors. As an additional side benefit, the simplicity of Eq. (44) enables us to extract closed-form expressions for the normal velocity component and its peak value. As outlined in Eqs. (19) and (20), continuity may be used to retrieve

As usual, the maximum value of $v$ occurs at the edge of the layer, namely,

Despite its simplicity, the ability of Eq. (45) to predict its Blasius analogue is quite satisfactory. Recalling from Eq. (22) that the normal velocity corresponding to the classic Blasius solution is proportional to *αa*/2, we get

It is reassuring that the prediction based on the optimized quartic profile is off by less than 1.6%.

### D. Extension to higher-order polynomial approximations

Pursuant to this analysis, higher-order polynomials can be constructed by securing additional boundary conditions and properties associated with the Blasius problem. For example, quintic, sextic, and septic polynomials can be judiciously constructed to mirror the behavior of the Blasius solution at the two endpoints of the viscous layer. These can be arrived at using a similar optimization technique, namely, by solving for the optimal values of *s*^{*} and other coefficients that stand to minimize the total residual of several additional constraints. To sketch this process, we start with *N* = 5 in Eq. (35) and then select three foundational conditions that must be fixed in the optimization procedure in order to properly capture the character of the Blasius solution. We find it essential to set

These immutable equalities enable us to determine three of the unknown constants, *a*_{0} = *a*_{2} = 0 and *a*_{3} = 0.99 − *a*_{1} − *a*_{4} − *a*_{5}. As for the remaining parameters, *a*_{1}, *a*_{4}, and *a*_{5}, they can be evaluated in such a way to minimize the objective function $R$ defined in Eq. (42). Thus, by simultaneously solving the system of constraints,

we extract *a*_{1} = *s*^{*} ≈ 1.630 398, *a*_{4} ≈ − 2.041 276, and *a*_{5} ≈ 1.291 600. Our optimized quintic polynomial becomes

Note that the value of *s*^{*} here has essentially converged onto the exact Blasius trajectory. Moreover, despite its simple four-term composition, Eq. (51) proves to be quite accurate, as it exhibits an *L*_{2} error of 0.002. In addition to providing a simple expression that can be conveniently integrated and differentiated, it reproduces the fundamental boundary-layer characteristics, precise in at least four decimal places, with a maximum relative error of 0.0012% in the estimation of *b*. Its incorporation within the KP approach leads to the following refined values, accurate in all digits displayed:

and

Extending this procedure further to a sextic polynomial can be achieved with relative ease. Using the same fixed conditions as before, the additional degree of freedom that we gain with *N* = 6 can be invested in minimizing an additional constraint such as $F\u2032(1)=F\xaf\u2032\u22480.090\u2009357\u2009643$ in Eq. (39). Our objective function becomes

As before, we first extract *a*_{0} = *a*_{2} = 0 and *a*_{3} = 0.99 − *a*_{1} − *a*_{4} − *a*_{5} − *a*_{6}. The remaining four constants {*a*_{1}, *a*_{4}, *a*_{5}, *a*_{6}} can be optimally determined such that

This system can be solved straightforwardly. One gets

With an *L*_{2} error of 0.001 758, the minor improvement that accompanies the sextic approximation is reflected in a maximum error of 0.0019%, which is accrued in the estimation of *b*; all characteristic properties remain, within the number of digits displayed, no different from those reported in Eq. (52). Hence, the use of an optimized sixth-order vs a fifth-order polynomial does not seem to be warranted, as the overall benefit remains minimal, except in the ability to better represent the endpoint conditions.

A septic polynomial, however, has the advantage of satisfying eight conditions identically, and these can absorb the values of the Blasius velocity function and its first three derivatives at both ends of the domain. In other words, a septic polynomial can be constructed either through optimization or in a manner to reproduce $F\u2032(0)=s\xaf$ precisely, along with seven other Blasius constraints. After some effort, we find that an assortment of eight physical requirements leads to a reasonably accurate solution. For example, using *N* = 7 in Eq. (35), the following eight conditions can be imposed:

and

These lead to *a*_{0} = *a*_{2} = *a*_{3} = 0 and $a1=s\xaf$, with the remaining constants yielding

As such, backward substitution into Eq. (35) returns

Despite its simple five-term composition, Eq. (58) proves to be quite accurate, with an *L*_{2} error of 0.002, which is comparable to the optimized quintic and sextic solutions. Nonetheless, although it matches the Blasius solution at the endpoints precisely, it reproduces the boundary-layer characteristics with a maximum error of 0.36%. For example, it predicts the following parameters and their corresponding relative errors (parenthetically):

Even at this high order, the benefit of optimization is clear. Instead of imposing the endpoint constraints on the third derivatives in Eq. (56), one may exclude *F*‴(0) = 0 and *F*‴(1) ≈ 4.4777 in favor of converting an additional Blasius condition into an objective constraint in Eq. (60). The latter may be augmented by the $F\u2032\u2032(1)=F\xaf\u2032\u2032\u2248\u22120.708\u2009537\u200962$ target requirement, namely,

As before, a systematic solution leads to *a*_{0} = *a*_{2} = 0 and *a*_{3} = 0.99 − *a*_{1} − *a*_{4} − *a*_{5} − *a*_{6} − *a*_{7}. This is followed by minimization with respect to {*a*_{1}, *a*_{4}, *a*_{5}, *a*_{6}, *a*_{7}}, which enables us to retrieve

In this case, the extremum corresponds to a zero residual for which $a1=s\xaf$ identically. Surprisingly, however, the *L*_{2} error of 0.0019 shows only a modest decrease relative to the sextic approximation. Nonetheless, by exhibiting an overall error that is comparable to that of the quintic solution, the septic expression is capable of restoring the boundary-layer properties with a maximum 0.0012% discrepancy in *β*. Other advantages stand in its ability to satisfy the endpoint values more accurately and to capture the boundary-layer properties more precisely than its five-term septic analogue [Eq. (58)]. This may be attributed to the latter being prescribed by the entire set of conditions in Eq. (56), notwithstanding the values of *α* and *β*. The inclusion of *α* and *β* as target functions in an optimization procedure is therefore essential.

### E. Improved effectiveness in thermal analysis

Having illustrated the utility of a refined formulation in a viscous-flow application, we can now proceed by applying the KP integral approach in conjunction with an energy balance to retrieve valuable estimates of the problem’s heat transfer and thermal boundary-layer characteristics. A sketch of our domain is provided in Fig. 3, where a section of the plate is heated, for *x* ≥ *x*_{0}, to a constant temperature *T*_{$w$}. As usual, the edge temperature *T*_{e} corresponds to the freestream conditions. In this setting, a thermal boundary layer, *δ*_{T}, will begin to develop starting at *x*_{0}. As shown by White,^{9} the flat-plate configuration can be modeled using

where *q*_{$w$}, *c*_{p}, and *k* are used to designate the wall heat flux, specific heat, and thermal conductivity, respectively.

To make further headway, an assumed velocity profile must be selected, say from Eqs. (28)–(30), and then paired with a shape-preserving temperature approximation. To draw a meaningful comparison, two cases will be considered: the first will apply Pohlhausen’s quadratic polynomial (with the smallest *L*_{2} error in Table I) and the second will repeat the analysis using the optimized quartic polynomial in Eq. (44). In the interest of brevity, the thermal boundary layer will be permitted to start at the leading edge, with *δ* = *δ*_{T} = 0 at *x* = *x*_{0} = 0 in Fig. 3. Using a similar quadratic polynomial to describe the temperature distribution within the thermal layer,^{9} we can write

When this approximation is substituted into Eq. (62), integrated over *δ*_{T}, and rearranged, one is left with

where *ζ* = *δ*_{T}/*δ* captures the thermal-to-viscous ratio of boundary-layer thicknesses. At this juncture, we may substitute $\delta (x)=ax1/2/U/\nu $ from Eq. (11), with $a=30\u22485.477$ from Table I, to obtain

where *α*_{T} = *k*/(*ρc*_{p}) is the thermal diffusivity. An asymptotic solution to this polynomial, valid over a practical range of Prandtl numbers, may be readily extracted. One finds

As shown with a dashed line in Fig. 4, our *O*(*Pr*^{−4/3}) series expansion agrees well with the numerical solution (solid line) as well as the simplest, one-term, *Pr*^{−1/3} expression (dotted line). With *ζ* and *δ* in hand, we can proceed to construct the local Nusselt number using Eq. (64). We get

In this case, the estimated *Nu*_{x} coefficient overshoots the correct value for constant wall temperature analysis of a flat plate by 10%, the correct constant being 0.332.

The classical KP procedure just described can be repeated using Eq. (30). First, Eq. (63) may be replaced by the optimized quartic polynomial estimate for the temperature such as

Inserting this expression into the energy balance integral in Eq. (62), we arrive at

Further expansions, integrations, and simplifications leave us with

or

The resulting expression can be put in the form $\zeta 3\u2212\alpha 0\zeta 5+\beta 0\zeta 6=\gamma 03Pr\u22121$ with $\alpha 0=27224$, $\beta 0=11480$, and $\gamma 0=3794203$. Here too, a robust approximation can be found, namely, $\zeta =\gamma 0Pr\u22121/3+13\alpha 0\gamma 03Pr\u22121\u221213\beta 0\gamma 04Pr\u22124/3+O(Pr\u22125/3)$. More visibly, we get

As illustrated with a chained (dashed-dotted) line in Fig. 4, this three-term perturbation series of *O*(*Pr*^{−5/3}) exhibits a steadily diminishing error with successive increases in the Prandtl number. It can be relied upon so long as *Pr* > 0.0045, knowing that its error drops precipitously below 1% for *Pr* > 0.45. Another “fractional” approximation that outperforms the *O*(*Pr*^{−5/3}) expansion in the 0.004 < *Pr* < 5.359 62 range can also be achieved using

These solutions and their relative errors are showcased in Fig. 4 along with the one-term theoretical *Pr*^{−1/3}. Despite the transcendental character of these equations, it is reassuring that they all return a root of *ζ* = 1 for *Pr* = 1, in compliance with the Reynolds analogy. Moreover, their differences become graphically imperceptible for *Pr* > 0.3 in Fig. 4(a) as their relative errors in Fig. 4(b) drop to 4% or lower. These observations confirm the viability of using *ζ* ≈ *Pr*^{−1/3}, which entails a less than 3.3% error, for all gases and fluids with *Pr* > 0.7.

To complete our thermal analysis in the wake of *ζ* and *δ* and recalling that $a=\delta Rex/x=9450379\u22484.9934$ from Table I, the local Nusselt number can be deduced systematically by taking

and so,

The resulting coefficient is only 0.5% higher than the “exact” coefficient of 0.332. As such, the use of an optimally constructed quartic profile *in lieu* of Pohlhausen’s is greatly beneficial. Clearly, simple integral techniques that are based on rationally conceived velocity and temperature profiles can yield valuable approximations for either local or global friction and heat transfer coefficients, with degrees of accuracy that range from less than 1%, as shown in this study, to 10% or more, depending on the precision of the assumed profiles.

In leaving this section, we note that although the same procedure can be extended to include pressure gradient and compressibility effects, such developments fall outside the scope of this work. Our interest will shift, instead, to the presentation of more accurate and continuously analytic solutions that remain uniformly valid as *ξ* → ∞.

## IV. COMPACT, UNIFORMLY VALID SOLUTIONS OF THE BLASIUS EQUATION

### A. Benchmark for flat-plate flow analysis

Shortly after the unveiling of the foundational boundary-layer equations by Prandtl,^{1} Blasius^{2} may have produced one of the most transformational similarity solutions in fluid mechanics. Although his conversion begins with the reduced Navier–Stokes equations, the ensuing work has proven to be an invaluable tool for the practical treatment of viscous flows. Not only has it provided deeper insight into boundary-layer effects; it has truly ushered in a series of similarity solutions and an era of unprecedented growth in the field.

In short, by introducing a similarity variable of the form $\eta =yU/(2\nu x)$, Blasius manages to reduce Prandtl’s boundary-layer equations to a simple nonlinear differential equation for the flat-plate flow problem outlined in Fig. 1. This equation would later serve as a permanent benchmark for testing new methods of analysis and a lasting beacon in the quest for new similarity solutions. Note that the factor “2” under the radical is adopted here for convenience. Then, using the same notation as before, we take $\psi =2\nu Uxf(\eta )$ to be the parental stream function. The corresponding velocities become

where overdots denote differentiation with respect to *η*; the latter is used to avoid confusion with *ξ*. Except for $F\xaf(\xi )=u/U=f\u0307$, dotted derivatives with respect to *η* differ from those referred to $\xi =\eta /\eta 99%=\eta 2/a$. When these transformations are substituted into Prandtl’s equations, a widely celebrated third-order differential equation is recovered, namely,

As for the boundary conditions, the velocity adherence at the wall, *u*(*x*, 0) = $v$(*x*, 0) = 0, and the freestream-merge condition, *u*(*x*, ∞) = *U*, readily convert through Eq. (75) into

### B. Series expansions for the Blasius problem

So far, the Blasius equation has only yielded exact solutions that are cast in the form of infinite series. Examples include the Homotopy Analysis Method (HAM) series by Liao^{47–49} and Liao and Campo,^{64} as well as the Adomian decomposition of an exponentially transformed Blasius variable by Ebaid and Al-Armani.^{65} Although Blasius^{2} provides asymptotic approximations for small and large *η*, his series expansion exhibits a finite convergence radius of 0 ≤ *η* < 4.023 464 493 5 beyond which the solution diverges. The Blasius expansion also requires knowledge of the initial gradient of the axial velocity with respect to *η*, denoted here as $\sigma =f\u0308(0)\u22480.469\u2009599\u200999$, a value that is often termed “connection parameter” or “Blasius constant.” Note that this value shifts to $\sigma /2\u22480.332\u2009057\u2009336\u2009215\u2009196\u20093$ when the similarity variable and stream function are written slightly differently as $\eta =yU/(\nu x)$ and $\psi =\nu Uxf(\eta )$. The corresponding Blasius equation becomes $f...+12ff\u0308=0$, with the same three conditions as in Eq. (77). Moreover, the initial slope changes to $s\xaf=F\xaf\u2032(0)=\sigma a/2\u22481.630\u2009398\u2009038\u2009629\u2009397$, when written as the wall derivative of the velocity function $F\xaf(\xi )=u/U$ of Sec. III A (cf. Table I).

Naturally, in view of its relevance to a wide variety of problems, the asymptotic behavior of the Blasius solution has been the subject of numerous investigations. Some of these inquiries have devoted themselves to its existence and uniqueness characteristics,^{66–69} including its tripole singularity inside its semi-infinite domain.^{70–72} Others have attempted to overcome its deceptive singularity using Padé approximants^{46,73,74} or Crocco variables that help relocate the singularity to a more convenient outpost.^{67,68,75}

In practice, the persistent efforts to solve the Blasius equation have led to vastly dissimilar algebraic expressions. These extend from simple piecewise approximations valid up to the edge of the boundary layer, i.e., 0 ≤ *η* ≤ *η*_{δ}, where $\eta \delta =\delta 99%U/(2\nu x)$,^{4} to infinite HAM series that remain valid over the entire range of 0 ≤ *η* < ∞.^{76} They have also given rise to a multitude of elegant approximations, such as those devised by Bairstow,^{44} Parlange, Braddock, and Sander,^{45} Boyd,^{50} and Iacono and Boyd.^{51}

Among those endeavors, several studies have devoted themselves to the determination of the Blasius constant with increasing degrees of success. For example, Bairstow^{44} and Goldstein^{77} reported 0.474 and 0.470, while Falkner^{78} and Howarth^{79} obtained 0.470 334 and 0.469 600. With the advent of modern computers, Fazio^{80} and Boyd^{72} arrived at 12 and 17 digits, which were later superseded by Abbasbandy and Bervillier,^{81} who achieved 21 digits of accuracy with their 0.469 599 988 361 013 304 509 value. The record for most significant digits at the time of this writing is held by Varin:^{82,83} he manages to secure 30 and then 100 digits while expressing $f\u0308(0)$ in the form of a convergent series of rational numbers to arbitrary order. In this study, the availability of the Blasius constant with a high degree of precision is taken into full consideration.

### C. Alignment of differential and integral predictions

In helping to establish the connection between the Blasius solution and the ongoing analysis of the KP approach, a well-resolved numerical solution of the Blasius equation is furnished in Table III. This is performed using the continuous Taylor series method to compute several characteristic parameters with sufficient accuracy in all digits displayed. Interestingly, the first numerical solutions of the Blasius equation were generated manually by Töepfer^{63} and, more precisely, by Howarth^{79} and then Cortell.^{84} Table III catalogs values of the characteristic stream function *f*, axial velocity $f\u0307$, and shear stress $f\u0308$. Other properties follow systematically. For example, since $f\u0307=0.99$ at *η*_{δ} ≈ 3.471 886 88, we are able to deduce the 99% boundary-layer thickness using

Note that the computed value for $a=2\eta \delta \u22484.91$ is about 1.8% lower than the traditional value of 5.0 computed manually^{63} and adopted still in most textbooks.^{9–12} As for the wall shear stress and momentum thickness, they can be readily evaluated from

and so,

Our refined Blasius calculations agree to varying degrees of accuracy with the coefficients associated with the integral approach; these have been specified earlier in Eqs. (11)–(14), particularly where they are estimated using one of the velocity profiles in Eqs. (28)–(30). Despite their basic agreement with the Blasius model, however, all piecewise approximations considered so far can be seen to suffer from two basic limitations. First, they all end abruptly at *y* = *δ*. Second, those considered up to the fourth polynomial order end with a value of *u* = *U* instead of *u* = 0.99*U*, thus slightly deviating from Prandtl’s 99% defining threshold at the boundary-layer edge. To overcome these inconsistencies, our attention is turned momentarily away from the analysis of piecewise approximations. Instead, the insight gained so far will be leveraged to develop uniformly valid continuously analytic profiles that mimic the behavior of the Blasius solution more expansively, i.e., in both the viscous region and the far field.

$\eta =yU2\nu x$ . | $\xi =\eta \eta \delta $ . | f(η)
. | $dfd\eta =F\xaf$ . | $d2fd\eta 2=f\u0308$ . | $d3fd\eta 3=f...$ . | $dF\xafd\xi =af\u03082$ . | $d2F\xafd\xi 2=af...2$ . |
---|---|---|---|---|---|---|---|

0.0 | 0.000 00 | 0.000 00 | 0.000 00 | 0.469 60 | −0.000 00 | 1.630 40 | −0.000 00 |

0.2 | 0.057 61 | 0.009 39 | 0.093 91 | 0.469 31 | −0.004 41 | 1.629 38 | −0.053 13 |

0.4 | 0.115 21 | 0.037 55 | 0.187 61 | 0.467 25 | −0.017 55 | 1.622 25 | −0.211 49 |

0.6 | 0.172 82 | 0.084 39 | 0.280 58 | 0.461 73 | −0.038 96 | 1.603 09 | −0.469 67 |

0.8 | 0.230 42 | 0.149 67 | 0.371 96 | 0.451 19 | −0.067 53 | 1.566 48 | −0.814 03 |

1.0 | 0.288 03 | 0.232 99 | 0.460 63 | 0.434 38 | −0.101 21 | 1.508 12 | −1.219 94 |

1.2 | 0.345 63 | 0.333 66 | 0.545 25 | 0.410 57 | −0.136 99 | 1.425 44 | −1.651 26 |

1.4 | 0.403 24 | 0.450 72 | 0.624 39 | 0.379 69 | −0.171 14 | 1.318 25 | −2.062 88 |

1.6 | 0.460 84 | 0.582 96 | 0.696 70 | 0.342 49 | −0.199 66 | 1.189 08 | −2.406 64 |

1.8 | 0.518 45 | 0.728 87 | 0.761 06 | 0.300 45 | −0.218 99 | 1.043 11 | −2.639 66 |

2.0 | 0.576 06 | 0.886 80 | 0.816 69 | 0.255 67 | −0.226 73 | 0.887 65 | −2.732 96 |

2.2 | 0.633 66 | 1.054 95 | 0.863 30 | 0.210 58 | −0.222 15 | 0.731 11 | −2.677 80 |

2.4 | 0.691 27 | 1.231 53 | 0.901 07 | 0.167 56 | −0.206 36 | 0.581 75 | −2.487 41 |

2.6 | 0.748 87 | 1.414 82 | 0.930 60 | 0.128 61 | −0.181 96 | 0.446 53 | −2.193 40 |

2.8 | 0.806 48 | 1.603 28 | 0.952 88 | 0.095 11 | −0.152 49 | 0.330 22 | −1.838 16 |

3.0 | 0.864 08 | 1.795 57 | 0.969 05 | 0.067 71 | −0.121 58 | 0.235 08 | −1.465 51 |

3.2 | 0.921 69 | 1.990 58 | 0.980 36 | 0.046 37 | −0.092 30 | 0.160 99 | −1.112 63 |

3.3 | 0.950 49 | 2.088 83 | 0.984 56 | 0.037 81 | −0.078 99 | 0.131 29 | −0.952 13 |

3.4 | 0.979 29 | 2.187 47 | 0.987 97 | 0.030 54 | −0.066 79 | 0.106 01 | −0.805 15 |

3.47† | 1.000 00 | 2.258 56 | 0.990 00 | 0.026 03 | −0.058 78 | 0.090 36 | −0.708 54 |

3.5 | 1.008 10 | 2.286 41 | 0.990 71 | 0.024 41 | −0.055 82 | 0.084 77 | −0.672 88 |

3.6 | 1.036 90 | 2.385 59 | 0.992 89 | 0.019 33 | −0.046 11 | 0.067 11 | −0.555 82 |

3.8 | 1.094 51 | 2.584 50 | 0.995 94 | 0.011 76 | −0.030 39 | 0.040 82 | −0.366 33 |

4.0 | 1.152 11 | 2.783 89 | 0.997 77 | 0.006 87 | −0.019 14 | 0.023 87 | −0.230 67 |

$\eta =yU2\nu x$ . | $\xi =\eta \eta \delta $ . | f(η)
. | $dfd\eta =F\xaf$ . | $d2fd\eta 2=f\u0308$ . | $d3fd\eta 3=f...$ . | $dF\xafd\xi =af\u03082$ . | $d2F\xafd\xi 2=af...2$ . |
---|---|---|---|---|---|---|---|

0.0 | 0.000 00 | 0.000 00 | 0.000 00 | 0.469 60 | −0.000 00 | 1.630 40 | −0.000 00 |

0.2 | 0.057 61 | 0.009 39 | 0.093 91 | 0.469 31 | −0.004 41 | 1.629 38 | −0.053 13 |

0.4 | 0.115 21 | 0.037 55 | 0.187 61 | 0.467 25 | −0.017 55 | 1.622 25 | −0.211 49 |

0.6 | 0.172 82 | 0.084 39 | 0.280 58 | 0.461 73 | −0.038 96 | 1.603 09 | −0.469 67 |

0.8 | 0.230 42 | 0.149 67 | 0.371 96 | 0.451 19 | −0.067 53 | 1.566 48 | −0.814 03 |

1.0 | 0.288 03 | 0.232 99 | 0.460 63 | 0.434 38 | −0.101 21 | 1.508 12 | −1.219 94 |

1.2 | 0.345 63 | 0.333 66 | 0.545 25 | 0.410 57 | −0.136 99 | 1.425 44 | −1.651 26 |

1.4 | 0.403 24 | 0.450 72 | 0.624 39 | 0.379 69 | −0.171 14 | 1.318 25 | −2.062 88 |

1.6 | 0.460 84 | 0.582 96 | 0.696 70 | 0.342 49 | −0.199 66 | 1.189 08 | −2.406 64 |

1.8 | 0.518 45 | 0.728 87 | 0.761 06 | 0.300 45 | −0.218 99 | 1.043 11 | −2.639 66 |

2.0 | 0.576 06 | 0.886 80 | 0.816 69 | 0.255 67 | −0.226 73 | 0.887 65 | −2.732 96 |

2.2 | 0.633 66 | 1.054 95 | 0.863 30 | 0.210 58 | −0.222 15 | 0.731 11 | −2.677 80 |

2.4 | 0.691 27 | 1.231 53 | 0.901 07 | 0.167 56 | −0.206 36 | 0.581 75 | −2.487 41 |

2.6 | 0.748 87 | 1.414 82 | 0.930 60 | 0.128 61 | −0.181 96 | 0.446 53 | −2.193 40 |

2.8 | 0.806 48 | 1.603 28 | 0.952 88 | 0.095 11 | −0.152 49 | 0.330 22 | −1.838 16 |

3.0 | 0.864 08 | 1.795 57 | 0.969 05 | 0.067 71 | −0.121 58 | 0.235 08 | −1.465 51 |

3.2 | 0.921 69 | 1.990 58 | 0.980 36 | 0.046 37 | −0.092 30 | 0.160 99 | −1.112 63 |

3.3 | 0.950 49 | 2.088 83 | 0.984 56 | 0.037 81 | −0.078 99 | 0.131 29 | −0.952 13 |

3.4 | 0.979 29 | 2.187 47 | 0.987 97 | 0.030 54 | −0.066 79 | 0.106 01 | −0.805 15 |

3.47† | 1.000 00 | 2.258 56 | 0.990 00 | 0.026 03 | −0.058 78 | 0.090 36 | −0.708 54 |

3.5 | 1.008 10 | 2.286 41 | 0.990 71 | 0.024 41 | −0.055 82 | 0.084 77 | −0.672 88 |

3.6 | 1.036 90 | 2.385 59 | 0.992 89 | 0.019 33 | −0.046 11 | 0.067 11 | −0.555 82 |

3.8 | 1.094 51 | 2.584 50 | 0.995 94 | 0.011 76 | −0.030 39 | 0.040 82 | −0.366 33 |

4.0 | 1.152 11 | 2.783 89 | 0.997 77 | 0.006 87 | −0.019 14 | 0.023 87 | −0.230 67 |

### D. Procedure to determine continuous velocity profiles

Based on several studies that seek to determine uniformly valid approximations of the Blasius equation,^{49,72,85,86} one may infer that an exponential basis function is appropriate to pursue. This choice is further guided by the slow decay of the Blasius solution in the far field. On this note, we proceed by adopting a function of the type, *F*(*ξ*) = 1 − exp[−*Y*(*ξ*)], where the structure of *Y*(*ξ*) → ∞ as *ξ* → ∞ is yet to be determined.

Then, based on our assessment of the piecewise solutions in Secs. III A–III D, we can infer that any refined formulation must seek to satisfy an amended form of Pohlhausen’s requirements in Eqs. (31)–(34), i.e., by imposing realistic conditions that do not restrain the solution unnecessarily or force it to stray away from the Blasius distribution.

This modification can be realized, first, by replacing the zero shear stress condition at *y* = *δ* by the actual velocity gradient at the wall, which plays a pivotal role in prescribing the model’s alignment with the Blasius form (being representative of the wall vorticity). Second, we insist on securing *u* = 0.99*U* at *y* = *δ*, in conformance with Prandtl’s strict definition of the boundary-layer thickness. Our amended conditions on *F* become

Note that we have eliminated Pohlhausen’s fifth condition, $F\u2032\u2032$(1) = 0, in favor of a fixed gradient at the wall that matches the Blasius constant, $s\xaf$. We have also omitted Pohlhausen’s fourth condition, which translates into *F*′(1) ≈ 0.0904. Although additional constraints can be adapted, we find this assortment to be sufficient at this order and consistent with the Blasius model requirements. In fact, by writing *Y*(*ξ*) as a cubic polynomial, we can systematically determine its coefficients based on Eq. (81). Using a polynomial basis for *Y*(*ξ*), we can put

A side benefit of this expansion is that the far-field condition, *F*(∞) → 1, becomes self-satisfied by the decaying exponential, granted that the highest-order coefficient *Y*_{N} is positive. The remaining effort is straightforward and rather interesting. First, the no-slip condition, *F*(0) = 0, returns 1 − exp(−*Y*_{0}) = 0 or *Y*_{0} = 0. Second, the fixed gradient, $F\u2032(0)=s\xaf$, yields $Y1=s\xaf$. Third, the vanishing momentum balance at the wall, $F\u2032\u2032$(0) = 0, requires $Y2=s\xaf2/2$. We are now one constant away from securing our objective. The last trailing constant may be determined from Prandtl’s formal definition of the boundary-layer cutoff point, namely, *F*(1) = 0.99. By writing $1\u2212exp[\u2212(s\xaf+12s\xaf2+Y3)]=0.99$, we retrieve $Y3=1.009\u200937s\xaf\u2248s\xaf$. Given the fortuitous nature of this outcome, we arrive at a surprisingly simple expression,

Evidently, with a cubic polynomial representation for *Y*(*ξ*), there is no freedom left to accommodate Pohlhausen’s adjusted fourth condition, *F*′(1) ≈ 0.0904, which corresponds to the asymptotically diminishing velocity gradient at the boundary-layer edge. Nonetheless, it is gratifying to see that Eq. (83) still agrees with the corresponding Blasius gradient within 3.2%.

### E. Comparison to other continuous solutions

In addition to Eq. (83) and the piecewise velocity approximations detailed in Sec. III A, several continuous models have been constructed for the purpose of capturing the Blasius characteristics. Of those, three analytic solutions for *u*/*U*, which continue to hold past the 99% disturbance thickness, will be considered and compared after expressing them in terms of *ξ*. These are

The foregoing approximations have been obtained using insight into the Blasius solution, rationalization, curve-fitting, or a combination thereof.

To start, the ability of Eq. (83) to reproduce the Blasius solution is confirmed in Table IV, where several values of *F* and its derivatives are evaluated and compared at both ends of the viscous layer. Compared to the other continuous profiles considered here, its predictions of *F*′(0), *F*(1), *F*′(1), and $F\u2032\u2032$(1) appear to be the closest to the Blasius values except in the case of $F\u2032\u2032$(1) ≈ −0.729, where the −0.722 curvature in the error function of Moeini and Chamani^{89} deviates from the Blasius value of −0.709 by 1.8%, whereas Eq. (83) differs by 2.9%. Moreover, its prediction of *F*‴(0) seems to be less precise than that of Savaş.^{88}

$uU=F(\xi )$ . | F(0)
. | F′(0)
. | $F\u2032\u2032$(0) . | F‴(0)
. | F(1)
. | F′(1)
. | $F\u2032\u2032$(1) . | F‴(1)
. |
---|---|---|---|---|---|---|---|---|

Blasius^{2} | 0 | 1.630 40 | 0 | 0 | 0.990 00 | 0.090 357 | −0.708 54 | 4.4777 |

$1\u2212e\u2212s\xaf\xi (1+12s\xaf\xi +\xi 2)$ | 0 | 1.630 40 | 0 | +1.1145 | 0.989 85 | 0.093 211 | −0.729 33 | 4.4753 |

$[tanhs\xaf\xi 5/3]3/5$ | 0 | 1.630 40 | 0 | 0 | 0.986 98 | 0.097 392 | −0.658 84 | 3.8398 |

erf(1.592 61ξ) | 0 | 1.797 07 | 0 | −9.1162 | 0.975 70 | 0.142 239 | −0.721 55 | 2.9397 |

$tanh(s\xaf\xi )$ | 0 | 1.630 40 | 0 | −8.6684 | 0.926 12 | 0.232 014 | −0.700 65 | 1.9404 |

1 − exp[−Y^{(4)}] [Eq. (86)] | 0 | 1.630 40 | 0 | 0 | 0.990 00 | 0.094 266 | −0.751 23 | 4.6264 |

1 − exp[−Y^{(5)}] [Eq. (88)] | 0 | 1.630 40 | 0 | 0 | 0.990 00 | 0.095 424 | −0.763 93 | 4.6676 |

1 − exp[−Y^{(6)}] [Eq. (89)] | 0 | 1.630 40 | 0 | 0 | 0.990 00 | 0.090 357 | −0.708 54 | 4.4678 |

1 − exp[−Y^{(7)}] [Eq. (90)] | 0 | 1.630 40 | 0 | 0 | 0.990 00 | 0.090 357 | −0.708 54 | 4.4777 |

$uU=F(\xi )$ . | F(0)
. | F′(0)
. | $F\u2032\u2032$(0) . | F‴(0)
. | F(1)
. | F′(1)
. | $F\u2032\u2032$(1) . | F‴(1)
. |
---|---|---|---|---|---|---|---|---|

Blasius^{2} | 0 | 1.630 40 | 0 | 0 | 0.990 00 | 0.090 357 | −0.708 54 | 4.4777 |

$1\u2212e\u2212s\xaf\xi (1+12s\xaf\xi +\xi 2)$ | 0 | 1.630 40 | 0 | +1.1145 | 0.989 85 | 0.093 211 | −0.729 33 | 4.4753 |

$[tanhs\xaf\xi 5/3]3/5$ | 0 | 1.630 40 | 0 | 0 | 0.986 98 | 0.097 392 | −0.658 84 | 3.8398 |

erf(1.592 61ξ) | 0 | 1.797 07 | 0 | −9.1162 | 0.975 70 | 0.142 239 | −0.721 55 | 2.9397 |

$tanh(s\xaf\xi )$ | 0 | 1.630 40 | 0 | −8.6684 | 0.926 12 | 0.232 014 | −0.700 65 | 1.9404 |

1 − exp[−Y^{(4)}] [Eq. (86)] | 0 | 1.630 40 | 0 | 0 | 0.990 00 | 0.094 266 | −0.751 23 | 4.6264 |

1 − exp[−Y^{(5)}] [Eq. (88)] | 0 | 1.630 40 | 0 | 0 | 0.990 00 | 0.095 424 | −0.763 93 | 4.6676 |

1 − exp[−Y^{(6)}] [Eq. (89)] | 0 | 1.630 40 | 0 | 0 | 0.990 00 | 0.090 357 | −0.708 54 | 4.4678 |

1 − exp[−Y^{(7)}] [Eq. (90)] | 0 | 1.630 40 | 0 | 0 | 0.990 00 | 0.090 357 | −0.708 54 | 4.4777 |

To be more specific, we are able to, first, corroborate the negative $F\u2032\u2032$(1) curvature of all velocity profiles, consistently with the −0.709 Blasius value at the edge of the viscous layer. The role of this negative curvature and its impact on Pohlhausen’s quartic polynomial are discussed extensively in Sec. III B. Among the uniformly valid profiles, the closest curvature corresponds to Yun’s hyperbolic tangent, while the farthest stems from Savaş’s hyperbolic tangent of fractional order. Second, we readily confirm that all models satisfy the no-slip and momentum balance requirements at the wall, where *F*(0) = $F\u2032\u2032$(0) = 0. Third, both the decaying exponential and Savaş profiles return *u*/*U* = 99% at *ξ* = 1, whereas Moeini–Chamani’s and Yun’s return 97.6% and 92.6%, respectively. As for the extraordinarily important velocity gradient at the wall, all three expressions by Savaş,^{88} Yun,^{87} and Eq. (83) restore the value of $F\u2032(0)=s\xaf\u22481.630$. Finally, as far as matching the shear value at the edge of the viscous layer, here too, both exponentially decaying Eq. (83) and Savaş profiles predict *F*′(1) ≈ 0.093 and 0.097, respectively, *in lieu* of 0.0904.

In short, Eq. (83) either matches identically the endpoint values of the Blasius distribution or deviates from them by less than 3.2% up to the second derivative in *F*. Naturally, an even higher-order approximation can be devised by taking *N* = 4, 5, … and securing additional conditions that are reflective of the Blasius solution, as performed in Sec. III D. When carefully constructed, all such solutions will have the intrinsic ability to approach unity as *ξ* → ∞, thus satisfying the far-field condition identically. They will also possess the freedom to accommodate additional endpoint requirements with successive increases in the polynomial order.

For further clarity, the four compact profiles in question are compared to the Blasius velocity function $F\xaf=f\u0307$ in Fig. 5(a) for 0 ≤ *ξ* ≤ 1.2 (0 ≤ *η* ≤ 4.1663), thus illustrating the degree by which they imitate the expected behavior both within the boundary layer and beyond. Graphically, it may be seen that the decaying exponential solution in Eq. (83) (dotted) is indiscernible from the Blasius curve. This may be viewed as being significant, given the relative simplicity of this profile. This is followed by the Savaş hyperbolic tangent of fractional order, which can be barely distinguished from the Blasius line, even in the magnified inset of Fig. 5(b). Conversely, Moeini–Chamani’s error function and Yun’s hyperbolic tangent are shown to deviate progressively, especially in the upper 0.5 ≤ *ξ* ≤ 1 portion of the viscous layer. In fact, a strong correlation may be seen to exist between the spatial agreement of a given flow approximation with the Blasius shape and its effectiveness at predicting the fundamental boundary-layer characteristics. This is confirmed in Table V, where boundary-layer predictions by the continuous profiles are cataloged and contrasted to the computed Blasius values. Unlike the results corresponding to the piecewise profiles in Table I, all integral properties that are posted here, such as *α* and *β*, are computed over the entire domain of validity, i.e., 0 ≤ *ξ* < ∞. As usual, the percentage error with respect to the Blasius benchmark, which accompanies each individual estimate, such as $\alpha ,\beta ,H,(\delta /x)Rex,CfRex$, $(\delta */x)Rex$, and $(v/U)\delta Rex$ is posted directly under it. For each profile, we also compute and display the overall $L2=\u222b01(F\u2212F\xaf)2d\xi 1/2$ error across the viscous domain to accurately quantify the cumulative discrepancy accrued in each continuous profile. Only the essential properties are tabulated here because other related quantities can be deduced fairly straightforwardly, e.g., $d=b=(\theta /x)Rex=CfRex$, $CDReL=2d$, and $G\u221e=[(\delta */x)Rex]/2$. The latter prescribes the maximum normal velocity, $(v/U)maxRex$, which occurs as *ξ* → ∞.

$F(\xi )=dfd\eta =uU$ . | $\alpha =\delta *\delta $ . | $\beta =\theta \delta $ . | $H=\alpha \beta $ . | $\delta xRex$ . | $CfRex$ . | $\delta *xRex$ . | $vU\delta Rex$ . | L_{2} error
. |
---|---|---|---|---|---|---|---|---|

Blasius | 0.350 47 | 0.135 26 | 2.591 10 | 4.909 99 | 0.664 11 | 1.720 79 | 0.833 40 | 0 |

$1\u2212e\u2212s\xaf\xi (1+12s\xaf\xi +\xi 2)$ | 0.351 | 0.136 | 2.586 | 4.904 | 0.665 | 1.720 | 0.833 | 6.4 × 10^{−4} |

Cubic model [Eq. (83)] (%) | 0.06 | 0.25 | 0.19 | 0.12 | 0.13 | 0.06 | 0.10 | |

$tanhs\xaf\xi 5/33/5$ | 0.351 | 0.136 | 2.574 | 4.891 | 0.667 | 1.716 | 0.824 | 2.0 × 10^{−3} |

Savaş^{88} (%) | 0.12 | 0.79 | 0.67 | 0.40 | 0.40 | 0.27 | 1.07 | |

erf(1.592 61ξ) | 0.354 | 0.147 | 2.414 | 4.949 | 0.726 | 1.753 | 0.800 | 1.6 × 10^{−2} |

Moeini–Chamani^{89} (%) | 1.08 | 8.49 | 6.83 | 0.80 | 9.35 | 1.89 | 4.02 | |

$tanh(s\xaf\xi )$ | 0.425 | 0.188 | 2.259 | 4.163 | 0.783 | 1.770 | 0.805 | 5.6 × 10^{−2} |

Yun^{87} (%) | 21.30 | 39.14 | 12.82 | 15.23 | 17.96 | 2.84 | 3.45 | |

1 − exp[−Y^{(4)}] | 0.352 | 0.136 | 2.587 | 4.896 | 0.666 | 1.723 | 0.836 | 1.7 × 10^{−3} |

Quartic model [Eq. (86)] (%) | 0.40 | 0.57 | 0.17 | 0.29 | 0.29 | 0.12 | 0.31 | |

1 − exp[−Y^{(5)}] | 0.352 | 0.136 | 2.586 | 4.892 | 0.667 | 1.723 | 0.837 | 2.2 × 10^{−3} |

Quintic model [Eq. (88)] (%) | 0.53 | 0.74 | 0.22 | 0.37 | 0.37 | 0.16 | 0.44 | |

1 − exp[−Y^{(6)}] | 0.350 | 0.135 | 2.591 | 4.911 | 0.664 | 1.720 | 0.832 | 2.5 × 10^{−4} |

Sextic model [Eq. (89)] (%) | 0.06 | 0.05 | 0.01 | 0.02 | 0.02 | 0.03 | 0.18 | |

1 − exp[−Y^{(7)}] | 0.350 | 0.135 | 2.591 | 4.910 | 0.664 | 1.721 | 0.832 | 1.1 × 10^{−4} |

Septic model [Eq. (90)] (%) | 0.02 | 0.02 | 0.01 | 0.01 | 0.01 | 0.02 | 0.15 |

$F(\xi )=dfd\eta =uU$ . | $\alpha =\delta *\delta $ . | $\beta =\theta \delta $ . | $H=\alpha \beta $ . | $\delta xRex$ . | $CfRex$ . | $\delta *xRex$ . | $vU\delta Rex$ . | L_{2} error
. |
---|---|---|---|---|---|---|---|---|

Blasius | 0.350 47 | 0.135 26 | 2.591 10 | 4.909 99 | 0.664 11 | 1.720 79 | 0.833 40 | 0 |

$1\u2212e\u2212s\xaf\xi (1+12s\xaf\xi +\xi 2)$ | 0.351 | 0.136 | 2.586 | 4.904 | 0.665 | 1.720 | 0.833 | 6.4 × 10^{−4} |

Cubic model [Eq. (83)] (%) | 0.06 | 0.25 | 0.19 | 0.12 | 0.13 | 0.06 | 0.10 | |

$tanhs\xaf\xi 5/33/5$ | 0.351 | 0.136 | 2.574 | 4.891 | 0.667 | 1.716 | 0.824 | 2.0 × 10^{−3} |

Savaş^{88} (%) | 0.12 | 0.79 | 0.67 | 0.40 | 0.40 | 0.27 | 1.07 | |

erf(1.592 61ξ) | 0.354 | 0.147 | 2.414 | 4.949 | 0.726 | 1.753 | 0.800 | 1.6 × 10^{−2} |

Moeini–Chamani^{89} (%) | 1.08 | 8.49 | 6.83 | 0.80 | 9.35 | 1.89 | 4.02 | |

$tanh(s\xaf\xi )$ | 0.425 | 0.188 | 2.259 | 4.163 | 0.783 | 1.770 | 0.805 | 5.6 × 10^{−2} |

Yun^{87} (%) | 21.30 | 39.14 | 12.82 | 15.23 | 17.96 | 2.84 | 3.45 | |

1 − exp[−Y^{(4)}] | 0.352 | 0.136 | 2.587 | 4.896 | 0.666 | 1.723 | 0.836 | 1.7 × 10^{−3} |

Quartic model [Eq. (86)] (%) | 0.40 | 0.57 | 0.17 | 0.29 | 0.29 | 0.12 | 0.31 | |

1 − exp[−Y^{(5)}] | 0.352 | 0.136 | 2.586 | 4.892 | 0.667 | 1.723 | 0.837 | 2.2 × 10^{−3} |

Quintic model [Eq. (88)] (%) | 0.53 | 0.74 | 0.22 | 0.37 | 0.37 | 0.16 | 0.44 | |

1 − exp[−Y^{(6)}] | 0.350 | 0.135 | 2.591 | 4.911 | 0.664 | 1.720 | 0.832 | 2.5 × 10^{−4} |

Sextic model [Eq. (89)] (%) | 0.06 | 0.05 | 0.01 | 0.02 | 0.02 | 0.03 | 0.18 | |

1 − exp[−Y^{(7)}] | 0.350 | 0.135 | 2.591 | 4.910 | 0.664 | 1.721 | 0.832 | 1.1 × 10^{−4} |

Septic model [Eq. (90)] (%) | 0.02 | 0.02 | 0.01 | 0.01 | 0.01 | 0.02 | 0.15 |

Starting with the *L*_{2} error, we compute 0.056 and 0.016 for the Yun and Moeini–Chamani profiles, respectively. These values are consistent with the level of spatial alignment between their curves and the Blasius distribution in Fig. 5. They are also reflected by their relative variations in estimating the boundary-layer properties, which range from 2.84% to 39.14% in Yun’s case and from 0.8% to 9.35% in Moeini–Chamani’s.

In hindsight, these levels are comparable to the effectiveness of the Pohlhausen polynomial profiles, with Moeini–Chamani’s results resembling those of Pohlhausen’s quadratic polynomial in Eq. (28). As for the overall *L*_{2} disparity, it falls to appreciably low levels of 2.0 × 10^{−3} and 6.4 × 10^{−4} for the Savaş profile and the decaying exponential relation in Eq. (83), respectively. These full orders of magnitude reductions may be viewed as being quite favorable because the corresponding percentage errors in estimating the principal boundary-layer properties suddenly decrease to virtually insignificant levels that do not exceed 0.79% for the Savaş profile and 0.25% for Eq. (83). These low levels can very well explain the reason for the decaying exponential form to become imperceptible from the Blasius solution in Fig. 5, including the graphically enlarged inset in Fig. 5(b). In fact, recalling that the Blasius equation itself is accompanied by a small truncation error that depends on the size of the flow Reynolds number and given that the Blasius model is only valid in the laminar range, the decaying exponential may be viewed as a simple, well-behaved, analytical solution of the Blasius equation in the 0 ≤ *ξ* < ∞ range. From an engineering perspective, however, the rationally optimized quartic polynomial, (5*ξ* − 3*ξ*^{3} + *ξ*^{4})/3, remains perhaps among the simplest polynomial profiles that can be conveniently integrated to produce closed-form expressions for the main boundary-layer properties. For the remaining class of continuous profiles, numerical evaluation is generally required to determine their integral properties.

### F. Extension to higher-order decaying exponential solutions

Using Eq. (82) and the cubic polynomial argument resulting in Eq. (83), it is possible to construct higher-order polynomial arguments that capture more precisely the conditions specified in Eq. (56). To illustrate this process, we begin with *N* = 4 in Eq. (82) and then apply five essential conditions, viz.,

and

After some algebra, a quartic polynomial for the exponential argument is obtained in the form of

where

Interestingly, due to the commonly shared physical requirements, the quadratic form of *X*^{(c)} and its corresponding cubic expansion *Y*^{(c)} will continue to recur at increasing orders.

In seeking a higher-order model with *N* = 5, the same procedure may be repeated with the addition of *F*‴(1) ≈ 4.477 704 1 to Eq. (85). The resulting quintic solution becomes

Note that only expansions with positive coefficients multiplying their highest-order terms are permitted to ensure that *Y*(∞) → ∞ and, as such, *F*(∞) → 1. For this reason, the construction of a sextic polynomial requires exchanging the *F*‴(1) condition with both *F*′(1) ≈ 0.090 357 643 and $F\u2032\u2032$(1) ≈ − 0.708 537 62, in addition to the five essential conditions in Eq. (85). The strict application of all seven constraints enables us to extract a sextic solution of the form

Finally, by imposing all eight conditions in Eq. (56), a highly accurate approximation can be achieved, namely,

The validity of these higher-order solutions is showcased in Table IV, where they are seen to mirror the behavior of the Blasius solution and its derivatives at the endpoints of the viscous domain. Therein, it may be ascertained that additional endpoint conditions are secured with each successive increase in the polynomial order, with the septic model satisfying all eight requirements identically.

The overall accuracy of these models is also affirmed in Table V, where their ability to recover the characteristic boundary-layer properties is systematically tested. In this vein, it can be seen that the relative errors associated with both sextic and septic solutions are very small, being *O*(10^{−4}), including their *L*_{2} values. For example, with an *L*_{2} error of 1.08 × 10^{−4}, the uniformly valid septic model appears to outperform most profiles in its class. As for the quartic and quintic models, their errors of *O*(10^{−3}) exceed those of their lower-order analogue in Eq. (83). In fact, it may be speculated that the additional nonlinearities that accompany both quartic and quintic models in Eqs. (86) and (88) do not warrant their further pursuit, as they do not seem to offer tangible advantages over their lower-order expansion. Conversely, the simple three-term decaying exponential in Eq. (83) stands as one of the most compact and accurate models in its class, being capable of surpassing its higher-order analogues in both overall precision and prediction of Blasius-related properties until the sixth order is achieved; the latter, however, is obtained at twice the algebraic cost and number of nonlinear terms. Note that none of the profiles beyond Eq. (83) appear in Fig. 5 as they become graphically indiscernible from the Blasius solution.

## V. CONCLUSIONS

In this study, we revisit a classic paradigm in fluid mechanics, namely, the Kármán–Pohlhausen (KP) momentum-integral approach, which has proven itself, time and time again, to be of tremendous academic and research value, especially in the analysis of viscous and thermal boundary layers. Our work is carried out under the auspices of four overarching themes that aim at achieving four seemingly distinct but effectively joint objectives.

The first aspiration, which has formed the initial motivation for this work, is the need to explain the paradoxical performance of Pohlhausen’s high-order polynomials when used in concert with the momentum-integral approach. As such, our first objective has focused on pinpointing the technical factors that have caused Pohlhausen’s quartic polynomial to accrue wider deviations during the evaluation of boundary-layer properties than its corresponding cubic and quadratic analogues. This perplexing behavior, we have found, is attributable to three specific boundary conditions, conceived by Pohlhausen in 1921, with the ability to compel the solution to stray from its expected outcomes.

Retrospectively, these conditions can be ranked in decreasing penalty levels. The first, which consists of a vanishing curvature at the boundary layer edge, $F\u2032\u2032$(1) = 0, is shown to be inconsistent with the asymptotic behavior of the exact solution due to Blasius. The latter predicts a negative curvature of $F\u2032\u2032$(1) ≈ − 0.709. The second, less conspicuous discrepancy is traceable to the vanishing shear requirement, *F*′(1) = 0, at the edge of the boundary layer, where the exact solution returns *F*′(1) ≈ 0.09. As for the third disparity, it is associated with the *F*(1) = 1 condition, which mildly disagrees with Prandtl’s strict definition of a boundary layer, i.e., as comprising the vertical distance to where *F*(1) = 0.99. Realizing that the last two requirements, no matter how different from their expected values, are actually observed by Pohlhausen’s quadratic and cubic polynomials, the reduced accuracy of the quartic profile relative to its lower-order expressions is pinned rather categorically on the prematurely imposed curvature requirement. Fortuitously, the systemic process of sorting and comparing Pohlhausen’s fundamental assumptions to the proper values of the Blasius velocity function and its derivatives across the viscous domain enables us to identify a more congruent set of boundary conditions. In fact, by imposing a slightly more realistic assortment of boundary conditions, another quartic expansion is arrived at, specifically *F* = (5*ξ* − 3*ξ*^{3} + *ξ*^{4})/3, which is accompanied by an order of magnitude reduction in error relative to the Blasius benchmark. This effort and its attendant optimization analysis evolve into a second objective; it also results in a systematic procedure for deriving rationally optimized velocity profiles at successive polynomial orders.

For example, while constructing an alternative quartic profile, we start by retaining four of Pohlhausen’s classic requirements; these include no slip and a judicious momentum balance at the wall, coupled with a smooth merging of the velocity and its gradient at the boundary-layer edge. We deem it essential to relax the zero curvature requirement and replace it with an objective function that accounts for the cumulative error residual affecting the most influential parameters in prescribing boundary-layer properties. Besides the velocity gradient at the wall, the two other parameters that we designate consist of the non-dimensional displacement and momentum thicknesses. By linking these two integral quantities to the polynomial coefficients representing the profile in question, the objective function at the fourth order is reduced to a sole function of the slope, *s*. The latter captures the wall vorticity and, as such, the initial velocity gradient. At length, by minimizing the residual with respect to *s*, the optimal initial trajectory *s*^{*} is found in a manner to produce the smallest variance from the Blasius results with a quartic profile.

The optimization procedure in Sec. III C, however, cannot be undertaken in isolation. To make further headway, linking the boundary-layer properties resulting from the KP integral approach to the actual polynomial coefficients in a velocity profile becomes, by necessity, our third objective. In the ongoing analysis, this requirement is met through the development of several reduced momentum-integral relations in Sec. II, and these are shown to greatly facilitate the evaluation of boundary-layer characteristics for a given velocity approximation. In fact, the resulting effort is found to be equally relevant to our first objective by virtue of its ability to enhance the effectiveness of tracing the source of deviation in any given flow profile, such as Pohlhausen’s quartic polynomial, to the functional form of the velocity under construction. In summary, several simple relations are provided at the conclusion of this effort, and these enable us to deduce the fundamental boundary-layer characteristics from three principal quantities: the initial velocity gradient at the wall and both displacement and momentum thicknesses in non-dimensional form, {*s*, *α*, *β*}. Properties that can be readily retrieved using Eqs. (11)–(22), (26), and (27) include the actual disturbance, displacement, momentum thicknesses, the skin friction and drag coefficients, the wall-normal velocity and, let us not forget, the vorticity. This explains why the effort to obtain the most precise set of {*s*, *α*, *β*}, from which all other properties can be derived, ends up dominating the optimization procedure.

In fact, as the optimization procedure is methodically extended to higher polynomial orders, we find it essential to retain only three physical requirements, particularly those that are matched identically by the Blasius model: *F*(0) = $F\u2032\u2032$(0) = 0 and *F*(1) = 0.99. The remaining conditions are added incrementally to our objective function with successive increases in the polynomial order and the corresponding number of undetermined coefficients. Logistically, by continuing to minimize the overall residual error, we actively ensure that the remaining boundary conditions, which are turned into objective constraints, are collectively secured, as closely as possible for the purpose of reducing deviations in boundary-layer predictions. This enables us to arrive in Sec. III D at the most effective quintic, sextic, and septic polynomials, which are capable of predicting the various boundary-layer metrics with an essentially 0.00% error relative to their Blasius values. We project, although we only show for the case of a quartic profile in Secs. III C–III E, that straightforward integration of an optimally constructed polynomial in viscous and thermal KP-type analyses will result in an overall improved framework. Naturally, since the error entailed in the KP approach is driven by the accuracy of the similarity-preserving velocity and temperature profiles that are adopted, the optimal refinement of the latter is due to favorably affect the former.

However, although polynomial approximations for velocity and temperature profiles are ideally suited for implementation in a global KP-type setting, where they facilitate the explicit evaluation of integral properties in closed form, they remain limited in their applicability to the boundary layer only. While comparing our optimal polynomials to a highly resolved Blasius solution, as necessitated by the above-stated objectives, we are reminded of their piecewise analytic forms, which end abruptly at *ξ* = 1. We also realize that extending their range of validity into the far field will allow them to mirror the Blasius distribution even as *ξ* → ∞, thus making them more realistic. Then, given the insight that we have gained into the merits of applying different endpoint requirements and recognizing that the Blasius far-field condition, *F*(∞) = 1, prohibits the use of a polynomial representation, we proceed to explore other options, such as the use of a decaying exponential for a more appropriate basis function, albeit at the expense of adding nonlinearities. The remaining endeavor, which focuses on pursuing uniformly valid solutions to the Blasius problem, becomes our fourth and final objective. Although many attempts have been made to solve the Blasius equation, most of them have been wholly numerical or resulting in infinite series solutions. Although there exist several integro-differential solutions that are formidably accurate,^{44–46,50,51} the complexity entailed in evaluating them prompts us to defer their analysis to future work.

As we redirect our investigation to the possibility of constructing uniformly valid solutions in Sec. IV, we exploit our understanding of the Blasius problem, which is developed while seeking to fulfill the other three overarching objectives, to identify five essential requirements that must be secured by a prospective formulation. These consist of *F*(0) = $F\u2032\u2032$(0) = 0, $F\u2032(0)=s\xaf$, *F*(1) = 0.99, and *F*(∞) = 1. The latter is identically satisfied by the basis function that we explore, namely, a decaying exponential, *F*(*ξ*) = 1 − exp[−*Y*(*ξ*)]. By casting *Y*(*ξ*) in the form of a positive polynomial in the far field, we are left with four constraints that can be applied to extract a surprisingly accurate and compact expression, $F(\xi )=1\u2212exp[\u2212s\xaf\xi (1+12s\xaf\xi +\xi 2)]$. By comparing this continuously analytic profile, which incurs an *L*_{2} error of 6.4 × 10^{−4}, to other models in its class, its viability is cemented as an equivalent solution to the Blasius problem in most practical applications. However, should greater precision be desired, higher-order polynomial arguments, including sextic and septic polynomials, that entail even smaller errors are also constructed. To the authors’ knowledge, the septic relation Eq. (90), which incurs an *L*_{2} error of 1.08 × 10^{−4}, may be representative of one of the most precise and uniformly valid analytical solutions of the Blasius problem using a relatively simple expression for the velocity itself and not one of its derivatives.

In closing, let us remark that this study, which was initially pursued in the context of streamlining the presentation of the KP momentum-integral analysis in a classroom setting, along with its connection to the Blasius problem, seems to have unraveled several useful findings. Given the appreciable academic interest in the KP approach, as well as the Blasius equation, we hope that the various explanations provided here, including the procedure to improve the accuracy of the KP method or to solve the Blasius problem using decaying exponential functions could be further explored.

With the advent of an optimal quartic solution that complements Pohlhausen’s cubic and quartic formulations, new challenges arise: the widely used *u* = *U*(2*ξ* − 2*ξ*^{3} + *ξ*^{4}) model has been incorporated into countless viscous and thermal studies that, it stands to reason, ought to be revisited using the more precise quartic model. This challenge can be daunting when taking into account that, in the absence of simpler alternatives over the course of a century, Pohlhausen’s quartic polynomial has become the staple of analytic approximations characterizing the KP approach. In a sense, it has formed the backbone of numerous laminar-flow solutions of viscous and thermal boundary layers, both with and without pressure gradients.^{27,28} Evidently, its reduced accuracy has prompted several researchers to seek better alternatives, such as the methods of Walz^{90} and Thwaites.^{91} However, knowing that the accuracy of the KP approach can be markedly improved when used in conjunction with more refined velocity approximations, it is hoped that its applicability may be extended to problems with variable pressure gradients and freestream velocities, where other methods have superseded its use. Other curve-fit techniques, such as the method of Thwaites, are successful at providing global predictions, albeit with no information about the actual velocity that is driving the solution from within the boundary layer. In contrast, the KP approach retains the advantage of associating any of its solutions with tangible velocity and temperature fields.

For this reason, we hope that this study, which increases our repertoire of boundary-layer approximations and solutions to the Blasius problem, will open up additional areas of research endeavor. Specifically, we hope that our findings will encourage future extensions of the KP approach, perhaps through the use of optimally constructed velocity approximations, to other compressible, viscous, and thermal flow applications.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## ACKNOWLEDGMENTS

This work was supported partly by the National Science Foundation through Grant No. 1761675 and the Hugh and Loeda Francis Chair of Excellence, Department of Aerospace Engineering, Auburn University. The authors are deeply indebted to Dr. Martin J. Chiaverini, Director of Propulsion Systems at Sierra Nevada Corporation, Madison, WI, for numerous technical exchanges and for his unwavering support of their flow field investigations. Finally, we profusely thank the Editor-in-Chief, Dr. A. Jeffrey Giacomin, and the anonymous reviewers for expediting the processing of this article.