A globally convergent and closed analytical solution of the Blasius equation with beneficial applications

The first similarity solution for the laminar boundary layer of a semi-infinite flat plate could be found in the work of Blasius in 1908,¹ i.e.

subjecting to Boundary Conditions (BCs)

And a series solution

was proposed, where

and

However, the series (5) converges only in a small domain [0, 1.8894/C), where C≈0.33206, therefore, the property of the flow near larger η could not be correctly mastered. Explicitly, the imperfection is caused by the trouble in implementing the boundary condition (4), which makes the solution (5) could not be closed at the infinite place. Probably attributing to this imperfection, some classical text books of fluid dynamics did not quote this solution.² 

Unfortunately, after more than a century, a globally convergent and Closed Analytical Solution (CAS) of the Eq. (1) is still beyond our scope, though many kinds of methods have been proved to be effective in approximately solving it.^3–18 Surely, for a practical purpose, sufficiently accurate Approximate Analytical Solutions (AAS)^{3–5,11–13,16,17} with relatively simple or compact forms could be satisfactory in using, however, they are after all not globally closed and could not completely reflect some overall features of the solution. Actually, the significance of finding a CAS of Blasius equation not only origins from ones’ desire for further understanding the laminar boundary layer but also from the necessities in treating other problems related with boundary layer, e.g. drag reduction,¹⁹ Sakiadis problems²⁰ or transversal flows coupled with laminar boundary layer flows,²¹ heat transfer^13,22 and boundary layer regulation.^23–26 

Therefore, even a small progress in this topic still seems meaningful. The efforts for solving the problem impelled progresses in numerical methods^{6–9,16,17,23} and approximate analytical methods.^10–15 Surely, more strength has been devoted in finding AAS of the equation, though finding a CAS still seemed impossible,¹⁸ because simple and compact analytical expressions looks more attractive. Perturbation methods¹⁰ or recently the homotopy analysis method^11–15 were both proved to be effective in finding AAS. Collocation methods with orthogonal function systems were also successfully tested in the problem.^16,17 Before the real CAS of the Blasius equation is found, these work provide powerful tools and effective insights. In this paper, therefore, we made our efforts for having a glimpse of the characteristics of its CAS.

After a series of trials, we guess that the difficulty for obtaining a CAS might not only attribute to the strong nonlinearity of the Eq. (1) and infinite definition domain but also to our overemphasizing the form of a CAS should own an explicit expression as f=f(η) or an implicit expression as F( f, η)=0, i.e. the relationship between f and η is directly established in one formula. Therefore, we guess the difficulty might be reduced if we abandon this burden but turn to finding an indirect expression as f=f(T), η=η(T), where T is a parametric variable. In this paper, the trial is performed and we find the solution of Eq. (1) indeed could be parametrically expressed by 2 series with very simple forms, whose convergence could also be easily determined. We organize the paper as follow. In section 2, by introducing 2 auxiliary variables, we reduce the Eq. (1) into a 2^nd order Ordinary Differential Equation (ODE). In section 3, we present a particular method for solving the order-reduced ODE and gives a theoretically closed and globally convergent analytical solution. In section 4, we will demonstrate some utilizations of the CAS and provide the analytical expressions of the thickness of the boundary layer, higher order derivatives of dimensionless velocity on the wall. The CAS was also utilized in analytically formulating the heat transfer in the boundary layer. The problem of boundary layer control or regulation through wall injection and slip velocity was also mathematically illustrated. Finally, the conclusion is given.

It is easy to transform the Eq. (1) to

therefore, we introduce the area coordinate as an auxiliary variable

and define

then we have

By the derivative transformation

Eq. (11) can be changed into

so that

and

According to the transformation

and Eqs. (10) and (11), then the Eq. (14) becomes

The form of Eq. (16) is actually equivalent to the equation derived based on Crocco transformation.²⁰ It is quite easy to prove that f, T and S are all positive and increasing in [0, +∞). Therefore, we can see in Eqs. (3), (4) and (9) that

Then the BCs (2) ∼ (4) of BE are transformed to the BCs of the ODE (16), as Table. I shows.

Based on the auxiliary variables, we have the useful transformations of the derivatives, e.g.

We propose a method for solving the Eq. (16) as follow. The method would also be applied to the problem of the boundary layer with wall injection and slip velocity.

Because R=1 at T=0, as Table I shows, we thus expand 1/R into a power series as

so that the Eq. (16) becomes

Easily, we can have the solution of the above equation

where {b_n} are coefficients of the power series (20) or (21), e₂ and e₃ are integral constants. According to Eq. (18), we have

In Table. I, we can see f=0 for T=0, therefore, according to Eq. (22), we have e₂=0. Similarly, we can see R=1 for T=0, so that e₃=1 according to Eq. (21).

Then the next work is to find the values of {b_n}. It is easy to see that the equality RR⁻¹=1 helps. Therefore, substituting series (20) and (21) into the equality gives

where

In Eq. (23), it is easy to see {b_n} are the functions of γ. By directly mathematical calculation, the first 30 items of {b_n} can be easily obtained, as listed in Table. II.

A closed solution must satisfy the final BC, i.e. the solution (21) must be closed at R=0 for T=1 (see Table. I). According to Eqs. (21) and (24), we then have

Theoretically, if we can obtain the real value of γ by solving the algebraic equation (25), then the solution of series (21) is completely closed. Unfortunately, there are infinite terms in Eq. (25), therefore, the real value of γ is impossible to be obtained. However, if the root of Eq. (25) ensures the convergence of the series (21), then the unknown function f and independent variable η could be both expressed by power series based on the relations (18) and (19): i.e.

It is easy to prove that Eq. (26) is a closed solution of BE (1). Given an arbitrary value of T, a correspondence between f and η is established. The following work is to prove that the solution (26) converges in the whole domain [0, 1).

By the mathematical induction, we can easily find the coefficients {b_n} obey a law as

where

The details of the operations could be found in the  Appendix A of this paper. We can approximately infer that the sequence {d_3k =a_3k/a_3k-3} might have a limit d. Please refer to the  Appendix B for the further discussion. However, the exact value of this limit seems impossible to be derived by pure mathematical analysis. Therefore, we turn to the numerical methods. By directly calculating the former 7590 items of a_n, we can have an insight into its distribution law, see Fig. 1. The calculation stops at k=2530, i.e. n_max=7590, because a₇₅₉₀=9.75×10³⁰⁷ nearly exceeds the maximum real number that a PC could express. Coupled with the mathematical analysis in  Appendix B, the calculation in Fig. 1 indeed indicates the existence of the limit of {d_3k} and the value of the limit is about

That is to say, for sufficiently large n=3k,

so that, for sufficiently large k, according to Eqs. (27) and (28), we have

In 0≤T<1, obviously, the function f and η are finite, i.e. the series (26) must converges to finite values. That is to say, for sufficiently large n, term b_nTⁿ<1, otherwise the convergence of the series (26) would be predominated by the harmonic series Σ1/n and diverges. Therefore we can reasonably infer that b_n≤1 for sufficiently large n, i.e.

Meanwhile, we know f and η are infinite at T=1, hence, series (26) must diverge at T=1. That is to say, for sufficiently large n, the term b_nTⁿ≥1 at T=1, or b_n≥1, otherwise the convergence of the series (26) would be predominated by the geometric series Σ(e^0.28γ)^k and converges. Therefore we can reasonably infer that b_n≥1 for sufficiently large n, i.e.

Synthesizing the both considerations, thus we are encouraged to infer that, for sufficiently large n,

Therefore, according to Eq. (29), we can approximately estimate e^0.28γ≈1, so that γ≈e^−0.28. And according to Eq. (24), we have

which is very close to the value 0.33206 obtained in prior studies.⁶ If we approximately solving the Eq. (25) by only keeping finite terms in the equation, i.e.

where ‘K’ indicates the order of approximation, then approximate values of C could also be obtained, as Table. III shows. The relative difference Er in Table. III is defined as

Obviously, though the value of C in Eq. (31) is derived by the convergence criterion (30), it is much accurate than the approximations in Table. III. Therefore, it is reasonable to infer that the convergence criterion (30) or (31) indeed approaches to the root of closure condition Eq. (25). In fact, an estimation based on pure mathematical analysis in  Appendix B also gave an accurate approximation about 0.3361 for C, which has a relative difference about 1.2%. Thus the reasonability of the analysis in this section could be relatively persuasive. In summary, a parametrical CAS as Eq. (26) indeed exists and converges globally, though an explicit expression directly connecting f and η is still in the mist. A direct calculation based on the solution (26) is provided in Fig. 2 with a comparison with the numerical result of L. Howarth.⁶ The accordance could indeed be found on the whole definition domain. Some tiny differences appear near T=1 in the Fig. 2b), because direct calculation is impossible to be performed there where infinite terms in Eq. (26) must be kept.

Comparing with approximate solutions, there are some natural advantages of the globally CAS (26), especially it could give some analytical results that have not been obtained in prior studies.

In this section, the higher order characteristics at the inner side (i.e. η=0) and outer side (i.e. η=∞) of the laminar boundary layer would be discussed. With the CAS (26), an accurate analytical expression for the higher order derivatives could be obtained on the wall or at T=0,

which prescribes the higher order characteristics of the curve of dimensionless velocity in Fig. 2b). Especially for sufficiently higher order derivatives, i.e. when n>>1, then

because of the convergence condition (30). Similarly, we can also have

In the Fig. 2b), the curve seems very flat near T=0 and the segment of the curve could even be regarded as straight. However Eq. (33) tells us that the property of the curve near T=0 deviates from our first impression. Meanwhile, Eq. (33) also shows that T=0 might be a higher order extremal point, which implies it might be a potential instability of the boundary layer.

When η approaches to the infinite or T approaches to 1, we can have an investigation on the difference η−f, whose physical meaning would be revealed directly related with displacement thickness of the boundary layer (see section 4.2). Employing the CAS (26), the difference is given

in which n=3k. Because the series Σ(n+1)⁻¹(n+2)⁻¹ converges and b_n→1 as given by the convergence condition Eq. (30), we then know this difference is limited. Explicitly, the limit in Eq. (36) of the difference η−f enables the BC (4) to be satisfied.

In former studies, boundary layer thickness δ* and δ** were given approximately or numerically.² Utilizing the CAS (26), analytical expressions of the boundary layer thickness could be given. Their definitions are given by

where y is the normal distance from the wall, which is given by

where

where u_∞ is the velocity of the incoming flow, ν is the kinematical viscosity and x is the longitudinal distance from the front tip of the flat plate. Then the integrals could be changed into

Employing the relation between η and T in Eq. (26), we can easily have

Substituting Eq. (38) into Eq. (37), the integrals could be calculated

Because the series Σ(n+1)⁻¹(n+2)⁻¹ converges, therefore analytical expressions of δ* and δ** in Eq. (39) are limited. In practice, expressions in Eq. (39) provide a method to calculate the thickness of boundary layer to an arbitrary accuracy. The convergences are also depicted in Fig. 3.

We can see the coefficients quickly approaches to values near 1.72 and 0.6638 respectively. Surely, these approximate values were ever obtained numerically or approximately before, however, this might be for the first time that they are analytically presented. Besides, we can further find the displacement thickness δ*, comparing with Eq. (36) and Eq. (39), is determined by the limit of η−f or Δ_∞ at the infinite:

Heat transfer in the boundary layer with a thermostatic wall is an important and familiar problem. Though seemingly simple the problem is, its analytical results are rare to be seen. The temperature is coupled with the flow in the boundary layer, which could be expressed as¹³ 

where Pr is the Prandtl number and θ is the dimensionless temperature of the fluid:

where H_w, H_∞ are temperatures of the wall and circumstance respectively and H is the temperature of the fluid in the boundary layer. In the past, the expression (40) could only be numerically calculated or analyzed by AAS¹³ even when Pr is a natural number. Fortunately, a theoretically feasible expanding could be performed by the CAS (26) if Pr is a natural number, because

such that

Because R is a polynomial with respective to T, as Eq. (21) shows, therefore Eq. (42) could be theoretically expanded and integrated analytically if the binomial theorem is employed. We are not going to do the tedious arithmetic operations for Pr>2, because it is purely arithmetic and is not essentially difficult in mathematics. Actually, three special situations, which could be analytically inspected individually, gain our special focuses.

The first situation is for Pr=1, then the expression Eq. (42) could easily give θ=T, i.e. dimensionless temperature θ changes in the same law as dimensionless velocity T in the laminar boundary layer.

The second situation is for Pr=2, whose analytical result is also rare to be seen in prior studies. However, based on the solution Eq. (21), the dimensionless temperature could be conveniently given by Eq. (42), i.e.

in which the denominator obviously converges. Simple numerical computation could indicate the denominator is almost totally dominated by the main part 1−γ/4. Therefore the expression of the dimensionless temperature in the boundary layer is expressed with respective to T.

Exceptionally importance belongs to the third situation, i.e. Pr=0, because it has been proved impossible to be analyzed even by AAS, e.g. method of homotopy analysis in the former studies,¹³ in which the situation of Pr=0 is indeed a singular point. The singularity makes the feature of the heat transfer for Pr=0 could only be obscurely glimpsed by the similar approximate analyses near Pr=0. Fortunately, the relation between Eqs. (20) and (21) provides the help, which naturally enables us to derive the exact results of Eq. (40) or Eq. (42) for Pr=0, i.e.

For T<1, the geometric sequence {Tⁿ} decreases much faster than the harmonic sequence {1/(n+1)}, therefore, the convergence of the series Σ Tⁿ/(n+1) is predominated by the geometric series Σ Tⁿ, namely, the numerator of Eq. (43) is finite for T<1. However, the denominator in Eq. (43) is infinite because b_n→1 as given in Eq. (30) and the harmonic series Σ(n+1)⁻¹ diverges, hence

Therefore, for 0≤T<1, θ=0 everywhere in the boundary layer for Pr=0. However, for T=1 or at η=∞, Eq. (43) gives θ=1, i.e. the fluid temperature equals the temperature of circumstance. The temperature gradient could also reveal the similar information. Without difficulties, a differential transformation could be given

Employing Eqs. (11) and (43) and because of Eq. (44), we can have

i.e. θ′=0 everywhere even at T=1 or η=∞. In summary, for Pr=0, we have

or

That is to say the temperature of the fluid H is homogeneous everywhere in the boundary layer and equals the temperature of the wall H_w, however the fluid at the infinite place still keeps the circumstantial temperature H_∞. None temperature gradient exists in the whole boundary layer.

The result in Eq. (48) seems odd and contradictory to our first impression, because a thermal boundary layer is usually regarded co-existing with the velocity boundary layer. However, this result actually accords well with the former observations. Though in former studies, a theoretical formulation of temperature distribution in boundary layer for Pr=0 was not given, the conclusions for 0<Pr<<1 has been given and confirmed that the heat transferring in the boundary layer is approaching to an independence from the flow.²² It thus can be seen the indications of Eq. (48) further confirms this assertion in a theoretical way. In fact, this is not difficult to be understood as long as we have a re-investigation on the physical meaning of Pr=0. We know the Prandtl number Pr is defined by

where μ, C_p and κ are the kinetic viscosity, specific heat at constant pressure and heat conductivity. Hence, Pr=0 implies the fluid has none heat capacity, i.e. C_p=0, or the fluid has infinite heat conductivity, i.e. κ=∞. This singular property makes the heat flux in the fluid instantly transfer to the infinite places with an infinitesimal temperature gradient. In another word, for this odd fluid, infinitesimal temperature gradient could establish an infinite heat flux because of the infinite heat conductivity κ and the heat convection is negligible comparing with the heat conduction.

These three situations has been plot in the Fig. 4. We are not going to the details for the analyses on other values of Pr, it is after all not the task of this paper. The work could be remained in the future, in which detailed expression for other Prandtl numbers could be provided and analyzed.

The significance of the work in this paper is not limited by presenting a CAS (26) for BE. Actually, the methods provided for solving the BE (1) and analyzing the convergence of the solution (26) supplied a feasibility to analytically investigate a class of problems related with the regulation of laminar boundary layer. By injecting fluid into or sucking fluid from the boundary layer^19,23–25 or introducing relative velocity between the fluid and flat plate with slippery surface,¹⁹ the boundary layer could be regulated. The problem is usually formulated as

The problem could be similarly processed as in section 2 for Blasius problem. We can obtain the same Equation as Eq. (16) but with BC correspondences different from Table. I. The new correspondences among the BCs are provided in Table. IV with the changes being shown in the cells with colored background.

Inhomogeneous BCs in Eq. (49) make the problem relatively more general and difficult than the Blasius Problem in Eqs. (1)∼(2). The augmented difficulty results in a rarity of analytical investigations on this problem.^19,21,23 Surprisingly, the method for solving Eq. (1), which was presented in the section 2 in this paper, could also be employed in the problem (49), although the mathematical complexity is significantly increased inevitably.

A new auxiliary variable should be introduced for transforming the interval [−σ, 1] of T into the interval [0, 1], because [0, 1] is more convenient for convergence analysis, therefore

in which τ belongs to [0, 1]. Consequently, the Eq. (16) is changed into

And the differential transformations in Eqs. (18) and (19) could be changed into:

Then the Eq. (51) could be solved as well as Eq. (16). For simplicity, denoting

and expanding 1/R into a series of τ:

then Eq. (51) becomes

i.e.

which gives

Utilizing the BC: R=1 for τ=0 in Table. IV, the integral constant w₄ is determined:

Because of Eq. (52) and Eq. (54), we have

then the integral constant w₃ could be determined by the BC: f=−ρ for τ=0, i.e.

The other coefficients {q_n, n≥1} of the series Eq. (55) could finally be determined by RR⁻¹=1. After an amount of tedious mathematical operations, we have

Finally, the solution Eq. (56) is closed by the BC at the infinite: R=0 for τ=1, i.e.

As can be seen in Eq. (54), we know q₀ is determined by C and σ. Therefore, given different values of σ and ρ, the friction feature C=f ″(0) of the boundary layer is regulated based on the relation Eq. (61). It is impossible to obtain the exact value of C from Eq. (61), because there are infinite terms in the equation. A rational choice is turning to solving this algebraic equation approximately by keeping sufficient amount of terms, i.e.

Based on the coefficients given in Eq. (60), approximate values of C for different given values of ρ and σ are determined as Fig. 5 shows, where N=1000. Because Eq. (62) is an approximation to the exact equation Eq. (61), therefore C≈0.328 for ρ=0 and σ=0, with a relative difference about 1.2% comparing with the value 0.33206. Dependences between the regulation parameters ρ and σ could be found on the contour lines of C. On the contour lines, ρ is generally decreasing progressively with respective to σ, which implies that the effect of wall injection ρ should be weakened as the effect of boundary slip σ is strengthened. The effects of larger values of regulation parameters ρ and σ were not contained in Fig. 5, because they introduce additional difficulty in solving the algebraic equation Eq. (62) with strong nonlinearity and this is not the aim of this paper, besides, they might break the prerequisite for establishing a steady laminar boundary layer of Blasius equation Eq. (1).

The work in this paper indicated that a series of nonlinear equations similar as Blasius equation (1) could be analytically solved and investigated. Although a Close Analytical Solution (CAS) must be parametrically expressed by two series, i.e. Eq. (26), it owns some natural advantages on analytically formulating higher order characteristics of boundary layer (see section 4.1) and the thickness of boundary layer (see section 4.2).

As a more general form of Blasius equation, Eq. (49) applies to more extensive problems related with boundary layers. From a pure mathematical perspective, Eq. (56) proposed a general solution to a class of nonlinear ODEs. Though we are also encouraged to give the expression of η, unfortunately, it is very complicated because of the existence of the denominator (1+σ)τ−σ in Eq. (53). Numerical integration might be simpler for calculating η here.

Some problem, which seems impossible to be analytically formulated in prior studies, e.g. heat transfer in the laminar boundary layer of a plate with constant temperature, could obtain some analytical expressions when Prandtl number Pr is a natural number. Especially for the singular situation with Pr=0, the CAS in this paper gives a satisfied answer. It is necessary to emphasize that other convection-diffusion problem coupled with laminar boundary layer, e.g. the concentration equation, could be similarly formulated.

A noticeable small mathematical trick should be reminded, i.e. expanding 1/R to a series first (see Eq. (20)) and next obtaining the series of R by Eq. (16), then determining their coefficients {b_n} by the relation RR⁻¹=1. However, if a usual method is adopted, e.g. firstly transforming the Eq. (16) into 2C²RR″+T=0 and only expanding R directly in it, then the subsequent convergence analysis would be very complicated, especially the heat transfer for Prandtl number Pr=0 (see Eq. (43)) would consequently be complicated, because the information between R and R⁻¹ was unintentionally ignored.

We thank the supports of National Natural and Science Foundation of China with the supporting number11502097, and Nature and Science Foundation of Jiangsu Province with the supporting number BK20130478, and Foundation of Senior Talent of Jiangsu University with the supporting number 1281130025. Gratitude is expressed to Dr. Chen Weiqi of the China Ship Scientific Research Center and Professor Yin Xieyuan of University of Science and Technology of China for their positive discussions. Many thanks for the beneficial advices of Dr. Song Mitao of Jiangsu University.

According to Eq. (23), explicitly, for k=1

And we could set

and assume

then substitute it into (23) we have

Because of Eq. (A1), then the expression in the right hand side of Eq. (A2) could be changed into

therefore

and

For the former K items in the sequence {a_n} or {a_3k}, i.e. a_3×1, a_3×2, a_3×3, ..., a_3K, two geometric sequences {c_3k} and {e_3k} could be found, which satisfies

where

where λ_c and λ_e are increasing rates. Then, we can have an estimation on a_3(k+1) by Eq. (A3), i.e.