On the K\'arm\'an momentum-integral approach and the Pohlhausen paradox

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(\xi)=u/U$. Here $\xi=y/\delta$ 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 behaviour, 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(\xi)=2\xi-2\xi^3 +\xi^4$. For example, our alternative profile, $F(\xi)=(5\xi-3\xi^3+\xi^4)/3$, leads to an order-of-magnitude improvement in precision when incorporated into the K\'arm\'an-Pohlhausen approach in both viscous and thermal analyses. Then noting the significance of the Blasius constant, $\bar{s} \approx 1.630398$, this approach is extended to construct a set of uniformly valid solutions, including $F(\xi)=1-\exp[-\bar{s}\xi(1+\bar{s}\xi/2+\xi^2)]$, which continues to hold beyond the boundary-layer edge as $y\rightarrow\infty$. Given its substantially reduced error, the latter is shown, through comparisons to other models, to be practically equivalent to the Blasius solution.


Introduction
Following Prandtl's pioneering treatment of boundary layers (Prandtl 1904), and the shapepreserving similarity solution reported by Blasius (1908), the Kármán-Pohlhausen (KP) momentum-integral approach, introduced through two sequential papers by Von Kármán (1921) and Pohlhausen (1921), is widely regarded as the staple in boundary-layer analysis. This is mainly due to its effectiveness at providing a wealth of useful detail 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, for example, those by Rosenhead (1963), Schlichting (1979), Cebeci & Cousteix (1998), Oleinik & Samokhin (1999), White (2006), Schetz & Bowersox (2011), Pritchard & Mitchell (2015), Schlichting & Gersten (2017), and others. It must be noted, however, that the KP approach is restricted to the analysis of the so-called inner region, with the outer far field referring to the 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 † Email address for correspondence: joe.majdalani@auburn.edu Waele power-law fluids (Bizzell & Slattery 1962), rarefied gases over porous plates (Jain & Kumar 1972), rotating three-dimensional boundary layers on cones (Bloor & Ingham 1977), thermal boundary layers on short rotating blades (Mohanty et al. 1977), forced convection from surfaces immersed in non-Newtonian fluids using planar ), axisymmetric , and irregular geometry (Nakayama & Koyama 1988), forced convection in highly pseudoplastic fluids (Andersson 1988), natural convection over bodies immersed in fluid-saturated porous media (Nakayama & Pop 1989), non-Newtonian motions within slider bearings (Bujurke & Jagadeeswar 1992), natural convection over bodies embedded in porous media (Mehta & Sood 1994), non-Newtonian boundary layers for low-Reynolds number mud flows (Huang & García 1998), radial film flows (Rao & Arakeri 1998), asymmetrically-mixed convective motions around horizontal cylinders (Mustapha & Mohand 2003), heat transfer from infinite isothermal cylinders with both Newtonian and non-Newtonian fluids (Khan et al. 2005(Khan et al. , 2006, momentum and thermal analyses of power-law fluids over flat plates (Myers 2010), and forced convection in power-law fluid-saturated porous media (Thayalan & Hung 2013).
In the present investigation, however, the object 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 its predictions and their true values for flat-plate motion. Our next objective is to revisit the assumptions made by Pohlhausen (1921), 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 (Pohlhausen 1921;Schlichting 1955), and those pursuing monotonically-increasing continuous profiles that merge smoothly into the far field as y → ∞ (Blasius 1908;Bairstow 1925;Parlange et al. 1981;Boyd 1997;Liao 1997Liao , 1999aBoyd 2008;Iacono & Boyd 2015). 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 labelled sketch of the streamline-bounded control volume and coordinate system in figure 1.
Since the majority of 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 which, when replaced by more suitable constraints, lead to substantially more accurate flow profiles. More specifically, we show that Pohlhausen's quartic polynomial, u/U = F(ξ) = 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 so doing, we arrive at an alternative representation, F(ξ) = (5ξ − 3ξ 3 + ξ 4 )/3, which will be shown to outperform other piecewise approximations in its class while predicting both viscous and thermal layer characteristics. It will also open several avenues of research inquiry, including higher-order extensions.
Recognizing the utility of this approach in deriving spatially-restricted piecewise approximations, and since no flat-plate study is complete without 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, or suction (see, for example, Levin et al. 2005;Zhao et al. 2017;de Paula et al. 2017;Ng et al. 2017;Hack & Moin 2017;Tsigklifis & Lucey 2017;Hewitt & Duck 2018;Tsang et al. 2018Tsang et al. , 2019Brynjell-Rahkola et al. 2019;Beneitez et al. 2019).
By taking advantage of the Blasius constant at the wall, F (0) =s ≈ 1.630398, 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 − exp[−sξ(1 + 1 2s ξ + ξ 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 all the way to infinity.
From an organizational standpoint, the paper is divided into three main sections, besides the introduction and conclusion. In § 2, 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 § 3 with a procedure to rapidly evaluate boundary-layer properties using piecewise velocity approximations. Along the way, the paradoxical behaviour 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 § 2 and 3 are extended in § 4 to the systematic development of continuous, decaying exponential profiles that mirror the characteristic behaviour 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.

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. These will play a key role in improving the effectiveness of the KP integral procedure, starting with the classical flat-plate flow problem that we revisit here as a test bed for analysis.

Flat-plate configuration
Our schematic of the physical setting in figure 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 labelled as δ. Streamlines situated outside this layer are deflected vertically by a distance δ*, which defines the displacement thickness. The upper streamline slings upwardly 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 uniform with u = U = const. It is useful to normalize the velocity and transverse coordinate using F(ξ) ≡ u/U and ξ ≡ y/δ.

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 = δ 0 (u/U) dy = δ 1 0 F(ξ) dξ. Then using ρ for density, a momentum balance enables us to write, as detailed by Schlichting (1979) or White (2006): where C f stands for the skin friction coefficient. As usual, 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. Next, we find it useful to define the non-dimensional displacement and momentum thicknesses, α and β, along with the shape factor, H, using whereα andβ represent the approximated integral values of α and β over the finite [0, 1] interval. Note the strict equalities ofα = α andβ = β for piecewise analytic profiles exhibiting the property F(ξ 1) = 1. For this reason, the distinction of usingα orβ 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. It may be also instructive to note that α, β and H are pure constants that depend strictly on the assumed velocity distribution F. † Since the focus of this study is on the steady form of (2.1), we continue with which, for a spatially-invariant freestream, becomes

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 (2.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 where µ 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 (2.6): These two expressions may be equated to obtain where ν is the kinematic viscosity. For a boundary layer that starts to develop at x = 0, one can use δ(0) = 0 and integrate (2.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 = U x/ν. 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 (2.6) leads to (2.12) † For uniformly valid, non-piecewise velocity profiles that continue to increase in the far field for ξ > 1, restoring the semi-infinite bounds on α and β in (2.2) and (2.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.
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 (2.12). Interestingly, we find the products and ratios of the traditional constants b and a, which are often tabulated in various monographs (Fox et al. 2004;White 2006;Pritchard & Mitchell 2015), to be proportional to s and β: Practically, for a given profile F, it is only necessary to calculate α, β and s to be able to deduce a = 2s/β, b = d = aβ and c = aα. As for the global drag coefficient, it may be evaluated using an elegant expression introduced by Von Kármán (1921), In view of (2.12), it is clear that C D = 2C 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 (2.11), the axial gradient of δ(x) can be eliminated in favour 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 To determine G δ explicitly, it is helpful to recreate the definition ofα in (2.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 δ = aα/2 = aα/2 = c/2 by way of (2.13). As such, and in reference to standard skin friction and momentum thickness tables (Schlichting 1979;Fox et al. 2004;White 2006;Pritchard & Mitchell 2015), the peak normal velocity coefficients, √ Re x (v/U) max , are simply half of the tabulated values for δ* √ Re x /x. However, for those approximations that observe Prandtl's F(1) = 0.99 cutoff requirement on u/U, we get G δ = aα/2 − 0.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 (2.20). Making use of F(∞) = 1, we can add and subtract the variable ξ by putting This simple outcome is, of course, contingent upon ξ F(ξ) − 1 → 0 as ξ → ∞. Such will be the case when The last condition in (2.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 (Fox et al. 2004;White 2006). 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 boundarylayer vorticity, ω. Being a cross-derivative function, the outcome is straightforward. We get where ω* denotes the non-dimensional vorticity. The last expression may be further simplified using (2.20) into ω* = −F (ξ) 1 + a 2 ξ 4Re x and so ω* max = ω* ξ=0 = −F (0) (2.27) Compared to the normal velocity and its peak value in (2.19) and (2.21), which are inversely proportional to the square root of the Reynolds number, the G -dependent term in (2.26) is appreciably smaller, being proportional to Re −1 x . 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 behaviour 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 (2.11) to (2.22), (2.26) and (2.27). In fact, the ability to connect various boundary-layer properties directly to an assumed velocity profile will prove essential in the remainder of this work.

Improved effectiveness of the KP 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 Von Kármán (1921), who provided the integral formulation, and Pohlhausen (1921), 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.

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 which are set equal to unity over the semi-infinite domain, 1 < y < ∞. These correspond to Pohlhausen's quadratic, cubic, and quartic profiles (Pohlhausen 1921): Schlichting's sinusoidal profile (Schlichting 1955): and a rationally-optimized quartic profile: Using their normalized forms, these profiles are collected in table 1 and compared to the classic solution due to Blasius (1908). Note that relative percentage errors are provided below each estimate along with the overall L 2 = [ 1 0 (F −F) 2 dξ] 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 (Schlichting 1979;Fox et al. 2004;White 2006;Anderson 2017;Pritchard & Mitchell 2015). Properties that are not shown can be easily deduced because, for laminar flow, (θ/x) which is already tabulated, and C D √ Re L is twice this constant. Conversely, the maximum wall-normal velocity, √ Re x (v/U) max , is simply half of the entries for δ* √ Re x /x, as shown in § 2.3. 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 behaviour 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 1, 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%.

Paradoxical behaviour of Pohlhausen's polynomials
So one may wonder, what is controlling the accuracy of these profiles? For several decades, 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 (1921) proposed five cardinal requirements which, apart from the wall adherence and smooth merging with the freestream conditions, i.e., the three constraints satisfied by the quadratic polynomial in (3.1), piled on two further requirements on the shear stress: one at the wall and one at the edge of the viscous domain. To summarize, the first three foundational constraints consist of: u(x, δ) = U and ∂u ∂y y=δ = 0 (smooth merging with the far field) (3.5) 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 µ ∂ 2 u ∂y 2 y=0 = dp dx (momentum balance at y = 0); dp dx = −ρU dU dx (Euler's far field) (3.6) Blasius (1908) 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 so do Schlichting's sinusoidal and the optimized quartic profiles. Lastly, 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 (Fox et al. 2004;White 2006). 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 boundarylayer coordinate ξ returns a curvature of −0.7085, which is not yet zero, but rather of order unity (figure 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 1). It also leads to a higher initial slope at y = ξ = 0, as one may readily infer from figure 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 figure 2(b). Interestingly, the assumption that u = U (and therefore F = 1 instead of 0.99) at y = δ (far right), can be seen to affect all of the piecewise approximations in figure 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 § 4.3.
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 2. 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 to 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 (3.7). In fact, it is only after recognizing these various limitations that one can manage to derive a rationallyoptimized 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 = δ.

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 behaviour is manifestly established in table 1, 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 figure 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 behaviour is desirable.
To begin, we opt for a non-dimensional velocity profile that can be approximated by an Nthorder polynomial of the form F(ξ) = N n=1 a n ξ n (3.8) Based on (2.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 (3.4-3.7) include, but are not limited to, On the one hand, when Pohlhausen's physical requirements (3.4-3.7) are applied on F(ξ), they translate into On the other hand, a robust numerical solution of the Blasius equation predicts, as shown in table 2, It can thus be seen that the F (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 (1) = 0 requirement must be released. The degree of freedom thereby gained can be leveraged while taking into account the direct relations identified in § 2.3, 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 (3.11) except F (1) = 0. Despite the penalty that the F (1) = 0 assumption engenders, it is retained here in observance of the Pohlhausen 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 (2.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 (Töepfer 1912), the classical Blasius solution predictedᾱ ≈ 0.344,β ≈ 0.133 ands ≈ 1.660 (Schlichting 1979;Fox et al. 2004). 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 § 2.3 to be the most influential in prescribing the remaining boundary-layer properties. At this juncture, we may substitute (3.14) into (3.15) and express R solely in terms of s. This enables us to minimize R(s) analytically, namely, by taking, dR ds s=s* = 0 where R(s) = 0.0514s 4 − 0.211s 3 + 0.833s 2 − 1.985s + 1.574 (3.16) Thus, 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 ofs ≈ 1.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. Lastly, by inserting this trajectory back into (3.13) 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 (1) = 0 constraint. One may also confirm that this solution leads to reliable estimates for the boundary-layer thickness and skin friction, namely, (δ/x) √ Re x ≈ 4.993 and C f √ Re x ≈ 0.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 1, its L 2 error of 0.008, evaluated over its finite interval, is another favourable metric: it reflects an order of magnitude reduction relative to its predecessors. As an additional side benefit, the simplicity of (3.17) enables us to extract closed-form expressions for the normal velocity component and its peak value. As outlined in (2.19) and (2.20), continuity may be used to retrieve As usual, the maximum value of v occurs at the edge of the layer, namely, In relation to the general form prescribed by (2.20) and (2.21), we recover here Despite its simplicity, the ability of (3.18) to predict its Blasius analogue is quite satisfactory.
Recalling from (2.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%.

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 behaviour of the Blasius solution at the two endpoints of the viscous layer. These can be arrived at using a similar optimization technique to the one pursued to obtain a quartic profile, 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 (3.8) 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 (3.15). Thus, by simultaneously solving the system of constraints, dR da 1 = 0 dR da 4 = 0 and dR da 5 = 0 (3.23) we extract a 1 = s* ≈ 1.630398, a 4 ≈ −2.041276 and a 5 ≈ 1.291600. Our optimized quintic polynomial becomes F(ξ) ≈ 1.6304ξ + 0.1093ξ 3 − 2.0413ξ 4 + 1.2916ξ 5 (3.24) Note that the value of s* here has essentially converged onto the exact Blasius trajectory. Moreover, despite its simple four-term composition, (3.24) 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: 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 (1) =F ≈ 0.090357643 in (3.12). 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 F(ξ) ≈ 1.6304ξ + 0.4950ξ 3 − 3.8897ξ 4 + 3.9920ξ 5 − 1.2377ξ 6 (3.28) With an L 2 error of 0.001758, 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 (3.25). Hence, the use of an optimized sixth-order versus 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 (0) =s 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 (3.8), the following eight conditions can be imposed: (3.29) These lead to a 0 = a 2 = a 3 = 0, a 1 =s, with the remaining constants yielding a 4 = 0.99 −s − a 5 − a 6 − a 7 a 5 = 3s − 2a 6 − 3a 7 − 3.86964 a 6 = 9.1843 − 6s − 3a 7 and a 7 = 10s − 16.7331 (3.30) As such, backward substitution into (3.8) returns F(ξ) ≈sξ − 1.83095ξ 4 + 0.93048ξ 5 + 0.68916ξ 6 − 0.42908ξ 7 (3.31) Despite its simple five-term composition, (3.31) 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 ( 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 F(ξ) ≈ 1.6304ξ + 0.6840ξ 3 − 5.0756ξ 4 + 6.5605ξ 5 − 3.5729ξ 6 + 0.7636ξ 7 (3.34) In this case, the extremum corresponds to a zero residual for which a 1 =s identically. Surprisingly, however, the L 2 error of 0.0019 does not 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 β. And though its overall error slightly exceeds that of the sextic solution, its main advantage stands in its ability to satisfy the endpoint values more accurately. It also captures the boundary-layer properties more precisely than its five-term septic analogue, (3.31), mainly because the latter is prescribed by the entire set of conditions in (3.29), notwithstanding the values of α and β. The inclusion of α and β as target functions in an optimization procedure is therefore essential.

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 figure 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 (2006), the flat-plate configuration can be modelled 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 (3.1-3.3), 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 1), and the second will repeat the analysis using the optimized quartic polynomial in (3.17). 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 figure 3. Using a similar quadratic polynomial to describe the temperature distribution within the thermal layer (White 2006 When this approximation is substituted into (3.35), 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 δ(x) = ax 1/2 / √ U/ν from (2.11), with a = √ 30 ≈ 5.477 from table 1, 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 ζ = 0.9283Pr −1/3 + 0.05745Pr −2/3 + 0.01067Pr With ζ and δ in hand, we can proceed to construct the local Nusselt number using (3.37). 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 (3.3). First, (3.36) may be replaced by the optimized quartic polynomial estimate for the temperature such as Inserting this expression into the energy balance integral (3.35), we arrive at The resulting expression can be put in the form ζ 3 − α 0 ζ 5 + β 0 ζ 6 = γ 3 0 Pr −1 with α 0 = 27 224 , β 0 = 11 480 , and γ 0 = 3 379 420 . Here too, a robust approximation can be found, namely, ζ = γ 0 Pr −1/3 + 1 3 α 0 γ 3 0 Pr −1 − 1 3 β 0 γ 4 0 Pr −4/3 + O(Pr −5/3 ). More visibly, we get ζ = 0.96634Pr −1/3 + 0.0362564Pr −1 + 0.00666116Pr −4/3 + O(Pr −5/3 ) ≈ Pr −1/3 (3.44) As illustrated with a chained (dash-dot) line in figure 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.35962 range can also be achieved using ζ = 1620Pr 1/3 − 9752Pr − 385 675Pr 1/3 − 3360Pr − 154 1 3Pr 1/3 (3.45) These solutions and their relative errors are showcased in figure 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 figure 4(a) as their relative errors in figure 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 = δ √ Re x /x = 9450 379 ≈ 4.9934 from table 1, the local Nusselt number can be deduced systematically by taking 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 ξ → ∞.

Benchmark for flat-plate flow analysis
Shortly after the unveiling of the foundational boundary-layer equations by Prandtl (1904), Blasius (1908) 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 η = y √ U/(2νx), Blasius manages to reduce Prandtl's boundary-layer equations to a simple nonlinear differential equation for the flat-plate flow problem outlined in figure 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 ψ = √ 2νU x f (η) to be the parental stream function. The corresponding velocities become where overdots denote differentiation with respect to η, which is deliberately used to avoid confusion with ξ. Except forF(ξ) = u/U =ḟ , dotted derivatives with respect to η differ from those referred to ξ = η/η 99% = η √ 2/a. When these transformations are substituted into Prandtl's equations, a widely celebrated third-order differential equation is recovered, namely, ...

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 (1997Liao ( , 1999a and Liao & Campo (2002), as well as the Adomian decomposition of an exponentiallytransformed Blasius variable by Ebaid & Al-Armani (2013). Although Blasius (1908) provides asymptotic approximations for small and large η, his series expansion exhibits a finite convergence radius of 0 η < 4.0234644935, 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 σ =f (0) ≈ 0.46959999, † a value that is often termed 'connection parameter' or 'Blasius constant.' Naturally, in view of its relevance to a wide variety of problems, the asymptotic behaviour of the Blasius solution has been the subject of numerous investigations. Some of these inquiries have devoted themselves to its existence and uniqueness characteristics (Weyl 1942;Callegari & Friedman 1968;Callegari & Nachman 1978;Fang et al. 2008), including its tri-pole singularity inside its semi-infinite domain (Punnis 1956a,b;Boyd 1999). Others have attempted to overcome its deceptive singularity using Padé approximants (Boyd 1997;Ahmad & Al-Barakati 2009;Peker et al. 2011) or Crocco variables to relocate the singularity to a more convenient outpost (Callegari & Friedman 1968;Callegari & Nachman 1978;Wang 2004).
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 η δ = δ 99% √ U/(2νx) (Pohlhausen 1921), to infinite HAM series that remain valid over the entire range of 0 η < ∞ (Liao 2010). They have also given rise to a multitude of elegant approximations, such as those devised by Bairstow (1925), Parlange et al. (1981), Boyd (2008), Iacono & Boyd (2015), and other dedicated researchers.
Among those endeavours, several studies have devoted themselves to the determination of the Blasius constant with increasing degrees of success. For example, Bairstow (1925) and Goldstein (1930) reported 0.474 and0.470, while Falkner (1936) and Howarth (1938) obtained 0.470334 and 0.469600. With the advent of modern computers, Fazio (1992) and Boyd (1999) arrived at 12 and 17 digits, which were later superseded by Abbasbandy & Bervillier (2011), who achieved 21 digits of accuracy with their 0.469599988361013304509 value. The record for most significant digits at the time of this writing is held by Varin (2014Varin ( , 2018: he manages to secure 30 and then 100 digits while expressingf (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 while seeking higher-order solution refinements. † This value shifts to σ/ √ 2 ≈ 0.3320573362151963 when the similarity variable and stream function are written slightly differently as η = y √ U/(νx) and ψ = √ νU x f (η). The corresponding Blasius equation becomes ... f + 1 2 ff = 0, with the same three conditions as in (4.3). Moreover, the initial slope changes tō s =F (0) = σa/ √ 2 ≈ 1.630398038629397, when written as the wall derivative of the velocity function

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 3. 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 (1912) and, more precisely, by Howarth (1938) and then Cortell (2005). Note that the computed value for a = √ 2η δ ≈ 4.91 is about 1.8% lower than the traditional value of '5.0' computed manually (Töepfer 1912) and adopted still in various textbooks (Batchelor 2000 46959999 (4.5) 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 (2.11-2.14), particularly, where they are estimated using one of the velocity profiles in (3.1-3.3). Despite their basic agreement with the Blasius model, however, all piecewise approximations considered so far can be seen to suffer from two basic limitations. Firstly, they all end abruptly at y = δ. Secondly, those considered up to the fourth polynomial order end with a value of u = U instead of u = 0.99U, thus slightly deviating from Prandtl's 99% defining threshold at the boundary-layer edge. To overcome these inconsistencies, we will turn our attention momentarily away from the analysis of piecewise approximations. Instead, we will proceed by leveraging the insight gained thus far to develop uniformly valid, continuously analytic profiles that mimic the behaviour of the Blasius solution more robustly and expansively, i.e., both in the viscous region and the far field.

Procedure to determine continuous velocity profiles
Based on several studies that seek to determine uniformly valid approximations of the Blasius equation (Boyd 1982(Boyd , 1999Liao 1999b;Parand et al. 2009), 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, , where the structure of Y(ξ) → ∞ as ξ → ∞ is yet to be determined. Then based on our assessment of the piecewise solutions in § 3.1- § 3.4, we can infer that any refined formulation must seek to satisfy an amended form of Pohlhausen's requirements in (3.4-3.7), i.e., by imposing realistic conditions that do not unnecessarily restrain the solution 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.99U 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 (1) = 0, in favour of a fixed gradient at the wall that equates to the Blasius constant,s. We have also omitted Pohlhausen's fourth condition, which translates into F (1) ≈ 0.0904. Although additional constraints can be imposed, 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 seek its coefficients based on (4.7). Using a polynomial basis for Y(ξ), we can put Y n ξ n ; N = 3 (4.8) A side benefit of this expansion is that the far-field condition, F(∞) → 1, becomes selfsatisfied by the decaying exponential, granted that the highest-order coefficient Y N is positive. The remaining effort is straightforward and rather interesting. Firstly, the no-slip condition, F(0) = 0, returns 1 − exp(−Y 0 ) = 0, or Y 0 = 0. Secondly, the fixed gradient, F (0) =s, yields Y 1 =s. Thirdly, the vanishing momentum balance at the wall, F (0) = 0, requires Y 2 =s 2 /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 − exp[−(s + 1 2s 2 + Y 3 )] = 0.99, we retrieve Y 3 = 1.00937s ≈s. 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, on the asymptotically diminishing velocity gradient at the boundary-layer edge. Nonetheless, it is gratifying to see that (4.9) still agrees with the corresponding Blasius gradient within 3.2%.

Comparison to other continuous solutions
In addition to (4.9) and the piecewise velocity approximations detailed in § 3.1, several continuous models have been developed with the goal of capturing the Blasius characteristics explored here. 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: Yun (2010) tanh (sξ) 5/3 3/5 Savaş (2012) erf(1.59261ξ) Moeini & Chamani (2017) (4.10) The foregoing approximations have been obtained using insight into the Blasius solution, rationalization, curve-fitting, or a combination thereof.
To start, the ability of (4.9) to reproduce the Blasius solution is confirmed in table 4, 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 (1) appear to be the closest to the Blasius values except in the case of F (1) ≈ −0.729, where the −0.722 curvature in the error function of Moeini & Chamani (2017) deviates from the Blasius value of −0.709 by 1.8%, whereas (4.9) differs by 2.9%. Moreover, its prediction of F (0) seems to be less precise than that of Savaş (2012).
To be more specific, we are able to, firstly, corroborate the negative F (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 § 3.2. 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. Secondly, we readily confirm that all models satisfy the no-slip and momentum balance requirements at the wall, where F(0) = F (0) = 0. Thirdly, 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ş,Yun and (4.9) restore the value of F (0) =s ≈ 1.630. Finally, as far as matching the shear value (or velocity gradient) at the edge of the viscous layer, here too, both exponentially-decaying (4.9) and Savaş profiles predict, respectively, F (1) ≈ 0.093 and 0.097 in lieu of 0.0904.
In short, (4.9) 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 higherorder approximation can be devised by taking N = 4, 5, . . . , and securing additional conditions that are reflective of the Blasius solution, as performed in § 3.4. When carefully constructed, all such solutions will have the intrinsic ability to approach unity as ξ → ∞, thus satisfying the Blasius (   [ t a n h ( ) ] s ‫ؘ‬ t a n h ( ) s ‫ؘ‬ e r f ( 1 . 5 9 2 6 1 ) ‫ؘ‬ Figure 5. Comparison of four continuous approximations for the Blasius solution according to Yun (2010), Savaş (2012), Moeini & Chamani (2017), and the systematically derived, decaying exponential function (4.9), which is imperceptible from the Blasius line. Results are shown (a) across the boundary layer into the far field, and (b) inside a designated quadrant where individual deviations from the Blasius curve are exaggerated.
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 functionF =ḟ in figure 5(a) for 0 ξ 1.2 (0 η 4.1663), thus illustrating the degree by which they imitate the expected behaviour both within the boundary layer and beyond. Graphically, it may be seen that the decaying exponential solution (4.9) (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 figure 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 For each profile, we also compute and display the overall L 2 = [ 1 0 (F −F) 2 dξ] 1/2 error across the viscous domain, to accurately quantify the cumulative discrepancy accrued in each continuous profile. As before, only the essential properties are tabulated because other related quantities can be deduced fairly straightforwardly, e.g., √ Re x ]/2. The latter prescribes the maximum normal velocity, (v/U) max √ Re x , which occurs as ξ → ∞.
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 figure 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 (3.1). In contrast, the overall L 2 disparity 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 (4.9), respectively. These full orders of magnitude reductions may be viewed as being quite favourable 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 (4.9). These low levels can very well explain the reason why the decaying exponential form is imperceptible from the Blasius solution in figure 5, including the graphically enlarged inset in figure 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, wellbehaved, 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, albeit confined to the viscous region only. For the remaining class of continuous profiles, numerical evaluation is generally required to determine their integral properties.
4.6. Extension to higher-order decaying exponential solutions Using (4.8) and the cubic polynomial argument resulting in (4.9), it is possible to construct higher-order polynomial arguments that capture more precisely the conditions specified in (3.29).
In seeking a higher-order model with N = 5, the same procedure may be repeated with the addition of F (1) ≈ 4.4777041 to (4.11). The resulting quintic solution becomes Y (5) (ξ) ≈sξ X (c) + 0.0523044ξ 3 + 0.0709987ξ 4 ≈ Y (c) + 0.0852769ξ 4 + 0.115756ξ 5 (4.14) 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.090357643 and F (1) ≈ −0.70853762, in addition to the five essential conditions in (4.11). The strict application of all seven constraints enables us to extract a sextic solution of the form Y (6) (ξ) ≈sξ X (c) + 0.418226ξ 3 − 0.350098ξ 4 + 0.0551749ξ 5 ≈ Y (c) + 0.681875ξ 4 − 0.570799ξ 5 + 0.0899571ξ 6 (4.15) Finally, by imposing all eight conditions in (3.29), a highly accurate approximation can be achieved, namely, Y (7) (ξ) ≈sξ X (c) + 0.335592ξ 3 − 0.106739ξ 4 − 0.183641ξ 5 + 0.078091ξ 6 ≈ Y (c) + 0.547148ξ 4 − 0.174026ξ 5 − 0.299408ξ 6 + 0.127319ξ 7 (4.16) The validity of these higher-order solutions is showcased in table 4, where they are seen to mirror the behaviour 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 5, 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 of 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 (4.9). In fact, it may be speculated that the additional nonlinearities that accompany both quartic and quintic models in (4.12) and (4.14) 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 (4.9) 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, albeit at twice the algebraic cost and number of nonlinear terms. Note that none of the profiles beyond (4.9) appear in figure 5 as they become graphically indiscernible from the Blasius solution.

Conclusion
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 behaviour, 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 (1) = 0, is shown to be inconsistent with the asymptotic behaviour of the exact solution due to Blasius. The latter predicts a negative curvature of F (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 compliant 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, evolves into a second objective; it also results in a systematic procedure to derive 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 the 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. Subsequently, 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 § 3.3, 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 § 2 that 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 (2.11) to (2.22), (2.26) and (2.27) include the actual disturbance, displacement, and 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 further 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, i.e., F(0) = F (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 § 3.4 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 § 3.3 and 3.5, 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 favourably 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 viscous 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 endeavour, 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 problem, most of them have been wholly numerical, or resulting in infinite series solutions. Although there exist several integro-differential solutions that are formidably accurate (Bairstow 1925;Parlange et al. 1981;Boyd 1997Boyd , 2008Iacono & Boyd 2015), the complexity entailed in evaluating them prompts us to defer their analysis to future work.
As we redirect our investigation into the possibility of constructing uniformly valid solutions in § 4, we exploit our understanding of the Blasius problem, which is developed while seeking to fulfil the other three overarching objectives, to identify five essential requirements that must be secured by a prospective formulation. These consist of F(0) = F (0) = 0, F (0) =s, 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(ξ) = 1 − exp[−sξ(1 + 1 2s ξ + ξ 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 provided. To the authors' knowledge, the septic relation (4.16), 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 unravelled 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 which, 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 (Khan et al. 2005(Khan et al. , 2006. Evidently, its reduced accuracy has prompted several researchers to seek better alternatives, such as the methods of Walz (1941) and Thwaites (1949). 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 can 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 to 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 endeavour. Specifically, we hope that our findings will encourage future investigations into seeking extensions of the KP approach, perhaps through the use of optimally-constructed velocity approximations, to other compressible, viscous and thermal flow applications.