Analytic Approximate Solution of the Extended Blasius Equation with Temperature-Dependent Viscosity

An explicit approximate solution is obtained for the extended Blasius equation subject to its well-known classical boundary conditions, where the viscosity coefficient is assumed to be positive and temperature-dependent, which arises in several important boundary layer problems in fluid dynamics. This problem extends a previous problem by Cortell (Appl Math Comput 170:706–710, 2005) when the viscosity is constant, in which a numerical solution was obtained. A comparison with other numerical solutions demonstrates that our approximate solution shows an enhancement over some of the existing numerical techniques. Moreover, highly accurate estimates for the skin-friction were calculated and found to be in good agreement with the numerical values obtained by Howarth (Proc R Soc A: Math Phys Eng Sci 164(919):547–579, 1938), Töpfer (Z Math Phys 60:397–398, 1912), and Cortell [34] when the viscosity is equal to 1, and when it is equal to 2.


Introduction
The Blasius equation was proposed by Blasius [1] in 1908 as a third-order nonlinear two-point boundary value problem to model the behavior of a steady state flow of viscous fluid. This is a two-dimensional incompressible laminar boundary layer 1 3 equation about a uniform stream of a non-Newtonian fluid over a semi-infinite flat plate. The semi-infinite plate is vertical and is placed parallel to the uni-directional flow, so the gradients of velocity and pressure in the x− direction can be ignored. Moreover, we usually require the velocity of the flow to be zero due to the stationary plate. This is called: the no-slip boundary condition. The classical Blasius laminar problem [2] reads Here, f represents the dimensionless stream function and is a dimensionless similarity variable.
The Blasius laminar problem is among the most important achievements in fluid mechanics because it laid the foundations for many principles and concepts that were established afterward. Blasius obtained a power series approximate solution for the equation, but the series wasn't converging on the whole domain but rather within only a finite interval. For this reason, finding a solution to the Blasius equation which is valid in the whole domain was of utmost interest to researchers. There is no closed-form exact solution for the equation, but many attempts have been made to find analytical approximate solutions and numerical solutions . Almost all methods and techniques were applied to obtain either numerical or approximate analytical solutions to the equation, Runge-Kutta method [3,11], series approximation and Padé approximant [2,17], numerical integration [3], perturbation [4,8], variational iteration method [5,14], homotopy method [6,23], Adomian decomposition method [9,13], recursive relations [10], finite difference method [12], shooting method [15,16], method of gradient [18], collocation method [19,21], generalized iterative differential quadrature method [20], combined Laplace transform and homotopy perturbation methods [22], numerical transformation methods [24], Haar wavelet approximation and a collocation method [25], optimal auxiliary functions method [26], direct algorithm method [27], quartic B-spline method [28], reproducing kernel method [29], optimal auxiliary functions method [30].
One critical value that characterize the problem is f �� (0) since the Blasius equation Eq. (1.1) can only be integrated directly to find f ′ if f �� (0) is known, then the problem can be converted to an initial value problem. So the value of f �� (0) was especially important, and it was the primary goal of Blasius in his work, and he obtained an estimated value between 0.3315 and 0.3317. In 1912 Töpfer [31] calculated to be 0.33206 using the Runge-Kutta method, which was confirmed by Howarth [3]. More accurate values were found numerically for the coefficient. Bairstow [32] calculated it as 0.33205735, Parlange et al. [2] and Girgin [20] obtained the value 0.33205734, and Fazio [33] obtained 0.332057336215, and Asaithambi [10] computed it to be 0.3320573362. The remaining results in the aforementioned research didn't give highly accurate values of . Liao [7] was the first to evaluate analytically, and it was found to be 0.33206 after computing the 35th-order series approximation. However, the 35th order approximate solution contained a huge number of terms that makes the approximation impractical and almost impossible to deal with.
A more general form of the problem is to consider the following extended version of the original Blasius problem where the coefficient a refers to the viscosity. Cortell [34] was the first to study the extended Blasius problem (1.3) subject to (1.2) for 1 ≤ a ≤ 2. He obtained numerical and complete solutions of the problem by using a Runge-Kutta algorithm and he found f �� (0) = 0.46960 for a = 1 and f �� (0) = 0.33206 for a = 2 , which is roughly close to the numerical values obtained previously.
The present paper investigates the extended Blasius problem with where a( ) > 0 is now a continuous function represents a temperature-dependent viscosity. It is well-known that dynamic viscosity decreases with the increase of temperature for liquids and increases with the increase of temperature for gases ([35], pp 27) and [36]. In engineering applications, the viscosity highly depends on the temperature [37,38]. This type of viscosity is important to nanofluids studies and helps obtain accurate results in heat transfer applications [39]. The influence of temperature-dependent viscosity on the laminar mixed free forced convective, boundary layer flow over a vertical semi-infinite isothermal flat plate has been investigated by researchers [39][40][41][42][43][44][45][46][47]. It is worth noting that the coefficient f �� (0) defines the skin friction around the plate, and describes the frictional force (or shear stress) generated by the fluid flow. Reducing skin friction is of fundamental importance in studying the turbulent boundary layers of liquid flows [47], and reducing the turbulent skin-friction drag force in air flows is a major concern in aerodynamic vehicles since reducing skin friction for it can affect the vehicle's emissions and fuel consumption and helps improve the efficiency of aircraft to fly faster higher, faster, and heavier [48], so no wonder it recently became a very active topic of research. The existence and uniqueness of the solution of (1.4) subject to (1.2) have been proved in [49]. Also, the numerical solutions to problems related to Blasius equations using the Adomian decomposition method can be found in [50,51]. The authors in [52] demonstrated that the solution for an arbitrary value of a can be obtained from the classical Blasius equation with a = 1 after a transformation. It was shown that for arbitrary a, f �� (0) = √ a where = f �� (0) = 0.46960 at a = 1.
In this paper, we aim to provide an analytic approximate solution in explicit form for the (1.4) subject to (1.2). In Sect. 2, we establish the explicit approximate solution. In Sect. 3, we analyze the solution and compare it with other numerical solutions.

3 2 Approximate Solutions
Rewrite the given nonlinear differential Eq. (1.4) into an equivalent form as follows since functions such f �� ( ) = 0 cannot be considered because they don't satisfy the boundary conditions. Integrating now Eq. (2.1) from 0 to , we obtain where f � (0) = 0 is taken into account, and integration by parts with f (0) = 0 yields which is a nonlinear Volterra integral equation of the first kind that will be used in the sequel for calculating the analytical approximate solution by a direct method. The constant C 1 = a(0)f �� (0) can also be determined by using the third boundary condition f � (∞) = 1 and Eq. (2.3) to get

Fundamental Lemmas
We need the following fundamental lemmas.
Given this we obtain the inequality (2.7). ◻

Upper and Lower Bounds of the Solution
The

Analysis of Solutions
We conclude here based on the following analysis that the best approximation to the BVP (1.4) with (1.2) is given by the lower bound approximation in the whole region 0 ≤ < ∞ with the corresponding boundary conditions. Thus We are now in the position to explore some mathematical results and investigate the numerical treatment of the boundary value problem (1.4) with (1.2). The goal is to compare between the approximate solution and the numerical solution in order to numerically validate the analytic solution. This can be done using different numerical techniques. In [53,54], the authors implemented Matlab and Runge-Kutta technique with interplant extensions to validate analytic solutions to the onedimensional Eyring-Powell fluid flow and to a Darcy-Forchheimer fluid flow. For our numerical validation, we use the Maple software, which is a very useful tool and provides sophisticated numerical techniques. On the other hand, with this software, we establish a useful program with simple statements to solve the problems of boundary value problems (BVP). The type of problem (BVP) is automatically detected, and an applicable algorithm is used. In our case, the used technique is the middefer method which is a midpoint method that uses the enhancement schemes. For the enhancement schemes, Richardson extrapolation is generally faster, but deferred corrections use less memory on difficult problems. Furthermore, this technique is capable of handling harmless end-point singularities that the trapezoidal scheme cannot. In addition, the numerical technique offers a benefits strategy with the use of the continuation method, which modifies the coefficient of the secondorder derivative of the equation. In this way, the global error is well decreased by an appropriate choice of the number of the maxmesh [55].
In the following, we will explore the proposed technique with an interesting case, where the coefficient a( ) is a constant, then we will extend the investigation to more complicated problems, where we consider a( ) as a function. In Fig. 1  On the other hand, the appropriate value of the constant C 1 = a(0)f �� (0) can be determined numerically, by calculating the value of the second derivative of f ( ) at = 0 . Thus, in order to reduce the numerical volume, we will give a fitting expression for the constant C 1 in terms of the coefficient a( ) . Before proceeding with this point, we will present some results for different values of the coefficient a.
In addition, we present in Table 1 a comparison between the solution and its derivative with the numerical corresponding calculations for the case a = 2.
To explore the role of the coefficient a( ) in the considered boundary value problem, we present the variation of the numerical solution and the analytic form of the solution for the case a( ) = 1 . In Fig. 2, we present the variation of the boundary value solution for different values of the constant C 1 and the numerical solution. We can see the good agreement between the numerical and analytical  solutions for the value of the constant C 1 = a(0)f �� num (0) = 0.4696 . In addition, we present in Table 2 a comparison between the solution and its derivative with the numerical corresponding calculations.
In Table 3 we present the analytic solutions and their derivatives with the corresponding numerical solutions for different values of the coefficient a( ) Now, we return to the point that is concerned with the choice of the constant C 1 . As previously cited, we use the expression C 1 = a(0)f �� (0) and substitute it with the second derivative from the numerical solution. Here we give a fitting expression for the constant C 1 in terms of the coefficient a( ) that can be used for any value of a without the need for the numerical value of the second derivative of the solution at the initial value.   Other interesting and instructive cases that can be considered here are concerned with the coefficient a( ) , which can be considered as a function. As the fluid viscosity is dependent on temperature, shear rate, and time, then we can assume it as a function and choose some analytic forms, in which the coefficient is increasing and decreasing function. We consider the first case that is associated with the case when the viscosity increases with the increase of , (i.e. when temperature decreases as in liquids). In this case we assume a( ) = e b , with b > 0 where following the lemma 2, we get for the constant 0 < C 1 = a(0)f �� num (0) , which can be fitted in terms of the constant b as C 1 = a(0)f �� (0) ≈ −0.0192b 3 + 0.1625b 2 + 0.603b + 0.3397 . In Table 4, we present the numerical solutions with the lower and the derivatives in terms of the independent variable for different values of the constant b = 0.5, 1, 2, 5 . The lower bound agrees very well with the numerical solution when the constant C 1 approaches the limit a(0)f �� (0).
The last useful case that can be investigated is concerned with the function a( ) = 1 n +b , where b is a constant. This case is associated with the case when viscosity decreases with the increase of eta (i.e. when temperature decreases as in gases). We present in Table 5, the numerical solution with the lower bound and the derivatives in terms of the independent variable with b = 1 and for different values of the constant n = 1, 2 , where its maximum can be determined from lemma 2 as C 1 = a(0)f �� num (0) . The lower bound is in excellent agreement with the numerical solution when the value of the constant C 1 approaches this maximum value a(0)f �� num (0).        Table 5 Comparison with numerical values for the case a( ) = 1 n +b . f ( ), f � ( ), and f �� ( ) are calculated from the Eq. 3.3, with C 1 = a(0)f �� num (0) and b = 1