Solving the Nonlinear Boundary Layer Flow Equations with Pressure Gradient and Radiation

: The physical problem under consideration is the boundary layer problem of an incompressible, laminar ﬂow, taking place over a ﬂat plate in the presence of a pressure gradient and radiation. For the mathematical formulation of the problem, the partial differential equations of continuity, energy, and momentum are taken into consideration with the boundary layer simpliﬁcations. Using the dimensionless Falkner–Skan transformation, a nonlinear, nonhomogeneous, coupled system of partial differential equations (PDEs) is obtained, which is solved via the homotopy analysis method. The obtained analytical solution describes radiation and pressure gradient effects on the boundary layer ﬂow. These analytical results reveal that the adverse or favorable pressure gradient inﬂuences the dimensionless velocity and the dimensionless temperature of the boundary layer. An adverse pressure gradient causes signiﬁcant changes on the dimensionless wall shear parameter and the dimensionless wall heat-transfer parameter. Thermal radiation inﬂuences the thermal boundary layer. The analytical results are in very good agreement with the corresponding numerical ones obtained using a modiﬁcation of the Keller’s-box method.


Introduction
The mathematical description of the laminar, incompressible boundary layer flow remains a challenging problem with several aspects still unexplored. The first to describe the steady two-dimensional boundary layer flow over a semi-infinite flat plate was Blasius [1]. The adverse pressure gradient on the laminar boundary layer flow was numerically examined by Howarth [2]. Other studies on the effects of adverse pressure gradients, heat transfer, or thermal radiation on the fluid flow include, although not limited to, [3][4][5][6][7][8]. The effect of radiation in the presence of a magnetic field or with temperature dependent viscosity was studied in [9,10].
At high temperatures, thermal radiation has a significant effect on the flow field, which has important applications in many areas of engineering. This effect has been reported in several works where various flows in the presence of radiation have been studied, such as the free convective laminar flow in [11,12] and [13], the flow of an optically thin gray fluid in [14,15], and [10] or the flow past a moving vertical plate or a vertical oscillating plate in [16,17].
In [18], a computational model was presented for convective and radiative heat transfer of very high temperature gas cooled and gaseous core reactors, with results that are in very good agreement with experimentally based correlations. In [19,20], the interaction of flow with radiation was studied in hypersonic boundary layer flows. In [21,22], a comprehensive numerical analysis of heat transfer in a rectangular enclosure was performed, and several effects on the fluid flow and heat transfer have been extensively studied, such as the effects of the Rayleigh number, thermal conductivity ratio, as well as internal surface emissivity. These analyses revealed a significant cooling effect of the internal volume together with an increase in the thermal conductivity ratio [22]. In [23], it was shown that the thermal behavior of a compressible turbulent flow was influenced by radiation and that the radiative heat flux plays a role in the recirculating zone. More recent studies show the impact of thermal radiation and the viscous dissipation effect on the dynamic three-dimensional rotating flow of carbon nanotubes [24]. In [25], the slip effects on magnetohydrodynamic (MHD) axisymmetric buoyant nanofluid flow over a stretching sheet in the presence of radiation and a chemical reaction were studied. In [26], the mixed convection stagnation-point flow of a nanofluid past a permeable stretching or shrinking sheet was numerically examined with suction, radiation, and a heat source/sink.
In this study, the standard formulation of the boundary layer equations is used and analytically solved until the point of flow separation, for the steady state, two-dimensional, laminar, boundary layer flow of an incompressible fluid over a flat plate under radiation, and pressure gradient effects. The nondimensional form of the corresponding equations with pressure gradient and thermal radiation are obtained in Section 2, utilizing the Falkner-Skan transformation [8]. The resulting nondimensional equations form a nonlinear, coupled system of partial differential equations (PDEs), accompanied by adequate boundary conditions. Such kinds of systems have been studied in the literature using mostly numerical methods since finding an analytical solution is very difficult, if at all possible. However, instead of using a numerical method in order to solve this system of PDEs, one may choose an approximate technique instead. The main advantage when using an approximate technique instead of a numerical is that a closed form approximate formula for the solution is obtained. This may lead to useful formulae, even if they are approximations, for quantities of physical interest. In addition, the implementation of approximate techniques is highly facilitated nowadays by the use of symbolic computation software, such as Mathematica or Maple.
Approximate methods have been successfully applied in the past in order to study fluid mechanics problems. The most well known approximate methods are the perturbation techniques [27,28]. For recent applications of perturbation techniques to fluid mechanics problems, one may refer to [29][30][31][32][33][34]. Other approximate techniques have also been applied in the study of various physical problems, including boundary layer flows. The references of the present paper are confined to few indicative papers of the last decade. The variational iteration method was utilized in [35] for solving the Blasius equation and in [36] for studying the magnetohydrodynamic flow over a nonlinear stretching sheet. The Adomian decomposition method was used not only for the solution of the Blasius equation in [37] but also for the study of micropolar flow in a porous channel [38]. Also, the differential transform method was used in [39] and [40] for the study of specific wedge flow and magnetohydrodynamic flow problems. It is worth mentioning that the differential transform method, combined with the method of steps, has been successfully applied for solving partial differential equations [41], as well as delayed differential equations and systems [42,43].
Maybe a breakthrough for the approximate methods was Liao's homotopy analysis method (HAM), which was introduced in [44]. Since then, it has been extensively used for a variety of physical problems, such as beam problems, oscillation problems, wave propagation problems, and of course, fluid mechanics problems. A thorough description of HAM can be found in the books [45,46], where it is also shown how other approximate methods can be unified by HAM. Probably the first papers where HAM was used for the study of boundary layer problems were [47][48][49], where the Blasius and the Falkner Skan equations were solved. Ever since, HAM has been successfully used, literally in hundreds of papers, for the study of boundary layer flows. Indicatively, it is mentioned that recently HAM has been used for the study of flows of nanofluids [24,[50][51][52][53][54][55] and flows in porous media [56][57][58][59].
In the vast majority of the papers regarding boundary layer flows where HAM was used, the nonlinear equations under consideration were homogeneous with boundary conditions at ∞. There are, however, also papers of this kind where the solved nonlinear equations are nonhomogeneous, but in most cases, the nonhomogeneous term is a constant. Even less are the papers of this kind, where the boundary conditions are not placed at ∞ but at some finite point.
In this paper, the steady state, two-dimensional, laminar, boundary layer flow of an incompressible fluid over a flat plate under radiation and pressure gradient effects, is solved for the first time using HAM, an analytical technique. This problem is mathematically described by a nonlinear, nonhomogeneous, and coupled system of PDEs. The nonhomogeneous terms in the momentum and energy equations are nonconstants, and in addition, the boundary conditions are placed not at ∞ but at the edge of the boundary layer. The appearance of such nonhomogeneous terms as well the specific type of boundary conditions, increase the difficulty of applying HAM. To the best of the authors' knowledge, only in [60] were similar equations and boundary conditions solved via HAM, where the two-dimensional, steady, boundary layer flow of a heat-conducting perfect gas was studied. However, no radiation effects were taken into consideration in the aforementioned paper.
The rest of the paper is organized as follows: In Section 2, the physical problem is given and the mathematical equations that describe it are derived. In Section 3, a short description of HAM is given together with its specific application for the physical problem studied in this paper. Finally, in Section 4, a thorough discussion of the obtained results is given with respect not only to the effect of a pressure gradient and radiation but also to numerical results concerning similar but not identical physical problems [3,5,6,19,20,61,62]. The presented results are novel, especially concerning the radiation effect on the boundary layer flow under the effect of a pressure gradient. More precisely, the obtained analytical results, reveal that the adverse and favorable pressure gradients substantially influence the dimensionless velocity of the boundary layer. Pressure gradient also plays a significant role in the changes of the dimensionless temperature of the boundary layer. Furthermore, the adverse pressure gradient causes significant changes on the dimensionless wall shear parameter and the dimensionless wall heat-transfer parameter, leading to boundary layer separation. Thermal radiation significantly influences the thermal boundary layer.

Mathematical Formulation
Consider the steady state, two-dimensional laminar flow of an incompressible fluid in the boundary layer that forms over a flat plate parallel to the free stream velocity, U ∞ . In a Cartesian coordinate system, Oxy, the surface is located at y = 0, 0 ≤ x ≤ L, as shown in Figure 1. The governing equations, after the boundary layer simplifications, form a parabolic system that constitutes of the continuity equation: ∂u ∂x the x-momentum equation: and the energy equation: with boundary conditions defined as follows: where u, v are the velocity components, T the temperature, p the pressure, ρ the density, ν the kinematic viscosity, µ the viscosity, C p the specific temperature of the fluid under constant pressure, k the coefficient of thermal conductivity, ∂q r ∂y the local radiant absorption, T w the temperature of the flat plate, u e , T e the velocity and the temperature at the edge of the boundary layer, respectively, and δ the boundary layer thickness.
Assuming that the temperature differences within the flow are sufficiently small, using the Taylor series of T 4 about T e and neglecting higher-order terms, T 4 can be expressed as a linear function of the temperature of the form . Then, the local radiant absorption is given by the expression: where α is the absorption coefficient and σ is the Stefan-Boltzmann constant [4]. At the edge of the boundary layer, Howarth's flow is considered, where the external velocity varies linearly with x, and at the edge, the velocity is u e = U ∞ 1 − x L , where U ∞ is the free stream velocity. By using Bernoulli's equation, the flow at the edge of the boundary layer and the pressure drop are associated: Moreover, introducing the dimensionless Falkner-Skan transformation and the dimensionless stream function [3]: as well as the dimensionless quantity the following system of PDEs is obtained: accompanied by the boundary conditions where the prime denotes partial differentiation with respect to η, Pr is the Prandtl number, f = f (x, η), θ = θ(x, η) andδ is again the boundary layer thickness, resulting now from Equation (7). Equations (8)-(10) are still dimensional. Thus by setting and dropping the · for reasons of simplicity, the following dimensionless equations are obtained 1 accompanied by the boundary conditions

A Short Description of HAM
As already mentioned, one of the most powerful approximate methods for solving nonlinear differential equations (DEs) is the Homotopy Analysis Method (HAM) as proposed by Liao, see [44][45][46]. According to HAM, instead of solving the original nonlinear DE, a series of linear DEs is solved and the solution of the nonlinear DE is expressed as a series of the solutions of these linear DEs. In the following, a short description of HAM is given, whereas in Section 3.2, it is demonstrated how it can be applied to System (12)- (14).
In order to keep things simple, HAM will be demonstrated by applying it to the nonlinear ordinary differential equation of n−th order accompanied by the initial conditions where y(t) is the unknown function and a i , i = 0, 1, . . . , n − 1 are known constants. According to HAM, an initial approximation y 0 (t) of y(t) should be constructed together with a homotopy H, so that • and y(t), i.e., the solution of Equations (15) and (16), is a solution of equation where q ∈ [0, 1] is a parameter, often referred to as embedding parameter. Liao proposed the homotopy [45]: where L is an auxiliary linear operator "inspired" by Equation (15) andh, H(t) are nonzero auxiliary quantites, called the convergence control parameter and auxiliary function, respectively. It is obvious that due to the linearity of L and which is Equation (15). This means that as q varies from 0 to 1, the solution Φ(t; q) of H[Φ(t; q), q] = 0, varies from y 0 (t) to y(t).
Since Φ(t; q) depends not only on t but also on q, it can be expressed in the form of a series as where satisfy the linear ODEs called high-order deformation equations, where and Thus, if all y k (t) can be found explicitly and the series on the right hand side of Equation (20) converges for q = 1, the solution of Equation (15) is given in the form of the so-called homotopy series In practice, K terms of the homotopy series are calculated, in order to obtain an adequate approximation of y(t). The corresponding series is called the approximate solution of K−order.
The HAM gives great freedom in choosing the initial approximation y 0 (t), the convergence control parameterh, the linear operator L, and the auxiliary function H(t). In order to choose these functions and operators, a set of base functions, for example is first chosen. This set S B is supposed to adequately describe the expected solution of Equations (15) and (16). Then, y 0 (t) should be chosen in such a way so that it a) can be expressed only in terms of S B and b) should satisfy the initial conditions (16), i.e., y 0 (0) = a 0 , y 0 (0) = a 1 , . . . y (n−1) 0 (0) = a n−1 .
The linear ODEs of Equation (21) are accompanied by the homogeneous initial conditions In this way, the homotopy Series (24) satisfies Equation (16). Then L and H(t) are chosen in such a way so that every y k (t), k = 1, 2, . . . is expressed as a series of e i (t), i.e., where c i,k , i = 0, 1, 2, . . . coefficients. In this way, y(t) is also expressed in the form where c i , i = 0, 1, 2, . . . coefficients. Moreover,h is chosen in such a way so that the series on the right hand side of Equation (24) is convergent. One way to chooseh, is to find the optimumh for which the so-called "discrete squared residual" error defined by is minimized, where t j = j∆t and ∆t is the step used to discretize the interval of t using N + 1. This has been proven to be a very good approximation of the total squared residual error; see, for example, [50] and [46] (Chapter 3). Another way to chooseh is to use the so-calledh−curves, which were proposed in [45] and have been used in a great amount of research papers. Theh−curves are graphs of quantities of y ap (t) calculated at specific values of t, which usually have a physical meaning. For example, y   (0) againsth, a horizontal line segment will be apparent, which is called the valid region ofh. If a value ofh is chosen from this region, it is quite certain that the Series (24) converges. Of course, maybe more than one important quantities with physical meaning can exist. Thus, more than oneh−curves can be used.

As initial approximations, the following functions
may be considered. It should be noted that the multiplicative term 1 1 − e −δ , as well as the term η in f 0 (x, η), are due to the boundary conditions atδ. The corresponding linear and nonlinear operators are chosen to be where the prime denotes partial differentiation with respect to η. Then, the higher-order deformation equations take the form and accompanied by the boundary conditions and where k = 1, 2, 3, . . ., χ k is defined by Equation (22), the auxiliary functions H f (η), H θ (η) are both chosen to be equal to 1, and the terms R f k and R θ k are determined via Equation (23) and are found to be: Thus, beginning with f 0 (x, η), θ 0 (x, η) and solving Equations (31), (33) and (32), (34), several approximations of f and θ can be found. It should be mentioned that the nonhomogeneous terms x u e du e dx and 16xασLT 3 e ρC p u e U ∞ , appearing in (12) and (13), respectively, result in the appearance of the terms k and R θ k , respectively. These terms greatly increase the complexity of Equations (31) and (32), and as a consequence, the computational time needed for obtaining several approximations of f and θ is also considerably increased. For example, in the case where both homogeneous terms for the pressure gradient and the radiation are considered, the computational time was 50% more compared to the homogeneous case, i.e., in the absence of the pressure gradient and/or radiation effect.
The proof of this theorem is omitted since it is similar to the proof of the general convergence theorem given in [45].
The values of the convergence control parametersh f andh θ are chosen in such a way so that the "discrete square residual" errors defined by are minimized for both f and θ, where η j = j∆η, ∆η = 0.16 and N = 50. Due to the fact that θ k does not appear in Equation (31), the optimum value forh f can be found independently and then used in order to find the optimum value forh θ . Table 1 shows the optimum values ofh f andh θ , as well as the corresponding errors for several values of the order K from 2 up to 12, in the case where L = 8, U ∞ = 1, x = 0.5, Pr = 0.7, T e = 300 K, a = 1/2, σ = 5.67 × 10 −8 , ρ = 1, c p = 1004 and |T w − T e | = 10 K. As the order increases, the error reduces for both equations, as expected. The values for the convergence control parameters are in accordance with theh-curves of the 15th order approximations, as shown in Figure 2. In this study, all the figures presented in the next section show results obtained utilizing the 15th order approximation for both dimensionless quantities, velocity, f , and temperature, θ. It is worth mentioning that in the case of no pressure gradient and/or no radiation effect, i.e., in the absence of nonhomogeneous terms, it suffices to consider only seven terms for f and only six terms for θ in order to obtain a good approximation of the solution.

Results and Discussion
To show the combined effect of a pressure gradient (adverse or favorable) and thermal radiation on the boundary layer flow, it is important to examine the dimensionless quantities of the problem under consideration, such as the dimensionless velocity, f (η), and the dimensionless temperature, θ(η). The shape of the boundary layer under adverse or favorable pressure gradient is clearly described by these quantities, and it turns out that it is very different from the shape of the boundary layer in the case with no pressure gradient.
In this study, the adverse pressure effect is achieved by changing the length, L (L = 5, 6, 8, 12) and retaining the free stream velocity, U ∞ constant. The free-stream temperature was considered to be T e = 300 K, Pr = 0.7, and c p = 1004. The temperature difference between the plate and the free stream was considered to be, |T w − T e | = 10 K. The radiation parameter, R, is given by the expression, where a is the absorption coefficient (a = 1/2), σ is the Stefan-Boltzman constant (σ = 5.67 × 10 −8 ), L is a characteristic length (the length of the flat plate), u e is the velocity at the edge of the boundary layer, U ∞ is the free stream velocity, c p is the specific heat under constant pressure, and ρ is the density of the fluid [10,61,62] (see Equation (13)). The pressure gradient effect plays a significant role on the flow, as indicated in Figure 3. The dimensionless velocity, f (η), is presented in Figure 3a. Adverse pressure gradient (APG) decreases the dimensionless flow velocity, while a favorable pressure gradient (FPG) increases the dimensionless velocity compared to the case with no pressure gradient (PG). The same effect is observed for the dimensionless temperature, θ (η), as shown in Figure 3b. The thermal boundary layer follows the behavior of the boundary layer. Numerically a similar problem has been studied without the effect of radiation [3,5,6].
A quantitative comparison of the analytical results obtained from HAM, with numerical results for the problem with the adverse pressure gradient, is presented in Figure 4. The numerical approach utilized is a modification of the Keller's-box method [5], where the governing equations are reduced to a first order system and then are approximated using central differences of second order accuracy. A multi-dimensional Newton-Raphson iteration scheme is then implemented for the solution of the resulting non-linear algebraic system, and the corresponding Jacobian matrix is not prescribed by the programmer but computed numerically. It can be observed that the numerical solution is in very good agreement with the analytical one, utilizing only 20 terms from HAM.  Additionally, the most important parameters for engineering applications are the skin friction coefficient, C f x and the local Stanton number, St x , defined per unit width of the plate, for the cases of a heating or cooling wall [5,6]. These parameters are functions of the dimensionless wall shear parameter, f w = f (x, 0), and the dimensionless wall heat-transfer parameter, θ w = θ (x, 0). Figure 5 shows the effect of adverse pressure gradient on the skin friction coefficient, f w with the dimensionless x-distance. It is observed that with different adverse pressure gradients, the dimensionless f w provides substantially different behavior. More precisely, as the APG is increased, the flow can separate faster concerning the x-direction, as indicated in Figure 5 for different realizations of the edge fluid velocity, u e = U ∞ (1 − x/L), representing the adverse pressure effect (ADV). Separation of the boundary layer, as shown in the ADV4 case, could lead to adverse effects and boundary layer separation.  Figure 6 shows the effect of APG on the dimensionless wall heat-transfer parameter, θ w with the dimensionless x-distance. It is observed that the dimensionless wall heat-transfer parameter θ w is appreciably influenced by the adverse pressure gradient, following the boundary layer behavior. This effect is more pronounced close to the separation region, as indicated in Figure 6. Finally, Figure 7 shows the effect of radiation on the thermal boundary layer. The flow is under the influence of an adverse pressure gradient. Close to the separation area (x = 0.8), a significant cooling effect of the thermal boundary layer is observed, which is mainly due to the influence of thermal radiation. This phenomenon was also observed in previous numerical and theoretical studies regarding compressible and turbulent flows [19,20,61,62]. Figure 7. Effect of radiation on the dimensionless temperature, θ (η), depending on the parameter, R, and the adverse pressure gradient, with respect to η.

Conclusions
The main goal of this work was to study the boundary layer problem of the incompressible, laminar flow past a flat plate in the presence of a pressure gradient and radiation. The physical problem is described by the continuity, the momentum, and the energy equations with the boundary layer simplifications. Under specific assumptions on the radiation effect and using the Bernoulli equation for the pressure drop and the Falkner-Skan transformation, a nonlinear, nonhomogeneous, and coupled set of PDEs is obtained, accompanied by boundary conditions at the edge of the boundary layer.
This system is solved using the HAM, although the specific form of the boundary conditions, increase the complexity of the initial approximations f 0 (x, η) and θ 0 (x, η), which together with the existence of the nonhomogeneous terms, increase the complexity of the higher-order deformation Equations (31) and (32). As a consequence, the computational time needed for obtaining several approximations of f and θ is considerably increased. It is indicative that in one of the cases, where both homogeneous terms for the pressure gradient and the radiation are considered, the computational time was 50% more compared to the homogeneous case. Moreover, the nonhomogeneous terms also affect the number of terms required in order to obtain a good approximation, which also increases the computational time. In some cases, only half terms are needed in the homogenous case.
To the best of the authors' knowledge, in only one previous study were similar equations and boundary conditions solved via HAM, highlighting the novelty of the current work. Additionally, no radiation effects were taken into consideration in the aforementioned paper. A short description of this mathematical method is given together with its specific application for the problem studied.
The obtained analytical results show that the pressure gradient plays a significant role in the changes of the dimensionless velocity and temperature of the boundary layer. An adverse pressure gradient causes significant alterations on important engineering parameters, such as the dimensionless wall shear parameter, f w , and the dimensionless wall heat-transfer parameter, θ w , leading to adverse effects and flow separation. An essential cooling of the boundary layer was found, revealing that thermal behavior is influenced by radiation. This finding is important because it could be utilized in many engineering applications, such as the high temperature, gas-cooled, and gaseous fuel nuclear reactors, as discussed in the introduction section. Thermal radiation could influence the flow field and thermal boundary layer rendering the above application as a possible control technique.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: