Received 20 April 2016; accepted 2 August 2016; published 5 August 2016
Ever since its first appearance in literature  the boundary layer theory under certain assumptions often leads to the application of similarity techniques. As a consequence, the independent spatial variables together with the unknowns can be converted from the original partial differential equations which are second order, nonlinear and coupled to an ordinary set provided that the flow domain and the boundary conditions are relatively simple.
The resulting set of ordinary differential equations rich in numerical and physical challenges and supplemented by boundary conditions can then be numerically addressed by means of initial value Ordinary Differential Equation (ODE) techniques. Two-dimensional flow over a surface creates a boundary layer as fluid particles move more slowly near the surface than near the free stream. With a similarity transformation, the boundary layer equation is converted to a set of ordinary differential equations known as the Blasius equation  or its more general counterpart, the Flakner-Skan (F-S) equation  . Since then, mathematicians and numerical analysts exposed to this huge challenge and motivated by a keen sense of exploration have delved into various facets of solving nonlinear boundary value problems through the determinism and simplicity of coupled one- dimensional initial value problems. The boundary layer equations adequately describe flows predicted by the F-S equations; however, the assumptions resulting from the similarity theory do not take care of the prescribed initial condition in general so the resulting solutions are accurate in some asymptotic sense (Ostrach  ). This explains why the value of the scalar must be adequately and iteratively accounted for as the independent variable approaches an infinite value. Continuous activity in this area has resulted in a large body of literature. Interested readers can assess a comprehensive list of these contributors in   .
The shooting algorithm is often regarded as a practical way to numerically handle the F-S boundary value problem. Many of the algorithms reported in literature require an initial guess of the shooting angle which should lead to acceptable results after an iterative numerical procedure. However an additional initial condition is required in order to replace the condition at infinity. To facilitate the numerical procedure, the initial boundary value problem is converted to a set of coupled nonlinear ODEs which are amenable to a straightforward iterative numerical solution by any of the ODE solvers like the Runge-Kutta method.
Boundary layer flow involving heat and energy transfer over a plate is very important in several engineering applications, for example, in polymer industry where plastic is produced or in general, for all cases involving the processing of sheet-like materials. In virtually all these cases, the sheet moves parallel to its own plane and induces the motion of the closest surrounding fluid as well as creates a boundary layer scenario. Some interesting physics involving the kinematics of the speed of the plate, stretching, contraction, heating and cooling have been observed to take place at this level (Shalini and Choudhary  , Bhattacharyya et al.  ). Earlier work involving flow past a stretching plate was recorded by Crane  . His work was further extended to deal with a permeable surface by Gupta and Gupta  . Cases involving unsteady shrinking film were considered by Wang  . This work was later vastly improved by Miklavčič and Wang  and Adhikary and Misra  . Similar work was carried out by Magyari and Keller  for self-similar boundary-layer flows induced by permeable stretching walls, as well as for stagnation point flow (Bhattacharyya  ). Shrinking sheet flow in this context, is regarded essentially as a backward flow which exhibits characteristics quite distinct from those of forward stretching flow.
This paper obtains the integral discrete analog of the energy equation by the way of the Green’s second identity and an appropriate complementary differential equation based on the Green’s function of the Laplace differential operator. Hybridization is achieved by retaining the integral discretization of the conventional boundary element method coupled with a finite-element type element-by-element determination of the scalar profile. This mechanism unlike the classical boundary element formalism permits “interaction” between the problem domain and its boundary by the intermediaries of its element based implementation as well as the accompanying interacting grid points. The resulting hybrid numerical procedure has therefore an implicit “local support” as well as an enhanced accuracy. These features are associated with both the Finite Element Method (FEM) and Boundary Element Method (BEM) by virtue of a straightforward numerical formulation as well as the subsequent element- driven implementation. We intentionally adopt this numerical route partly to address the scarcity of information on the boundary integral solutions for flows over a flat plate and also to seek for simpler numerical techniques to handle the species transport component of boundary layer flows.
2. Mathematical Formulation
The goal of the present investigation is to present a robust algorithm capable of producing faithful results for the numerical solution of the Falkner Skan (F-S) equation and accompanying energy equation. The application of boundary integral numerical techniques for the computation of F-S equation accompanied by heat and/or mass transfer is rare in scientific literature. Its emergence however is buoyed by the desire to seamlessly link the problem boundaries with its domain in order to relax the “boundary-only” approach of the BEM technique.
Consider two dimensional flow of an incompressible viscous fluid with heat transfer over a flat plate. The x-axis is taken along the plate while the y-axis is considered normal to the plate. In order to simplify the numerical solution, the following assumptions are made namely:
1) the fluid flow is considered laminar, stable and at steady state;
2) all body forces are neglected and all thermal properties are not temperature dependent.
The boundary layer equations can now be reduced to (Chow  ):
In the above system of equations, u and v represent velocities in the x and y directions respectively, T is the fluid temperature is the fluid density and the ratio is the kinematic viscosity where represent the fluid viscosity and p is the pressure. As long as the boundary layer is relatively thin, the external flow, and the pressure gradient will be independent of the thickness of the boundary layer. Equations (1) to (2) constitute the so called Blasius problem and it has been known that by stretching the vertical coordinate according to some law, the dimensionless velocity profiles measured in a laminar boundary layer at different distances x from the leading edge collapse into one. This paves the way for the so called self-similar solutions. The exter-
nal flow velocity and the pressure gradient are given by: where U is the exter-
nal velocity and b is a function of the flow geometry. The pressure gradient for a flat plate can be made favorable or adverse by assigning a positive or negative value to the variable m. Equations (1a)-(1c) are transformed by applying the non-dimensionalizations adopted by Falkner and Skan  which are similar but not identical to those of Blasius  . Non-dimensionalized flow coordinate and stream function and flow velocities yield:
A nondimensional stream function is found from the dimensional stream function
The velocity components become
A governing equation for f can be found by substituting these non-dimensionalized variables into the momentum equation to yield
Equation (6) is the well known Falkner-Skan equation. It comes with the following boundary conditions:
and reduces to the Blasius equation for. In the case of accelerating and decelerating flows respectively. Most physically realistic flows exist in the region of.
Over the years there have been several solutions of the F-S equation, it is not possible to go through all of them in depth. Among them are a host of numerical techniques involving finite differences, finite elements, spectral methods, adomian’s polynomials and perturbation methods. A comprehensive list of these can be found in Parand et al.  . In the calculations reported herein we adopt the shooting technique coupled with the Runge-Kutta ODE numerical solution method as the preferred numerical solution technique for the nonlinear ordinary differential equation.
We start by assuming a shooting angle say and iterate until the condition. The steps can be summarized as:
1) Guess initially and start evaluating in increments of the step size from zero to a.
2) If at a particular value of, then is decreased and is evaluated until for some. This mimics the bisection method and we can go ahead and bisect this range of values.
3) Once this range is identified, we continue with the iteration until oscillates around unity. We continue checking the trajectory and accuracy of our iteration comparing with a specified error tolerance.
4) We can then further refine our computation by applying the secant method.
Having then obtained the value of f and the derivatives as functions of we can go ahead to determine the required velocity profiles in the boundary layer.
Equation (1c) describes the energy aspect of the boundary layer flow. The assumption of incompressibility decouples the energy equation from the continuity and momentum equations.
The boundary conditions for velocity and temperature fields are given as
In order to solve the energy equation, the velocity components u and v are first obtained from the momentum and continuity equations and then substituted in terms of their similarity variables and f into the energy equation. We non-dimensionalize the energy equation by introducing a dimensionless temperature difference as:
where are the ambient temperature and the free stream velocity respectively, and is the specific heat at constant pressure. After substitution into Equation (1c) we obtain:
where is the Prandtl number and k is the thermal conductivity of the fluid. The transformed boundary conditions are:
Equation (10) paves the way for the introduction of a hybrid boundary-integral domain-discretized numerical procedure for a typical Blasius problem and its variations. Our major goal here is to achieve its integral representation by the way of the Green’s second identity.
The first step in this route is to propose a complementary differential equation of the type:
Equation (12) is interpreted physically to represent a one-dimensional diffusion subjected to a dirac-dellta type unit input. This characterizes a response known as the free-space Green function:
where K is an arbitrary constant, representing the longest element in the problem domain, and are field and source nodes respectively. The first derivative of G with respect to is
where is the Heaviside function defined as
The resulting integral equation yields:
As can be observed, Equation (16) is a typical Boundary Element Method (BEM) equation. However, its element-driven numerical implementation is the essence of the hybrid procedure adopted in this study. The problem domain to is divided into elements (Figure 1). Adopting a Finite Element Method (FEM) analog, Equation (16) can now be expressed as a summation of the integral representation of each element. For an M number of elements, Equation (16) can be written as:
Figure 1. Discretized problem domain.
Equation (18) is the interpretation of Equation (16) in elemental context where the variable “e” stands for an element and and are the coordinates of the end points of each element in the problem domain. We hasten to comment that in a typical element, there are two unknowns at each node namely the dependent variable and its spatial derivative. For the numerical implementation of Equation (18), we shall evaluate the source node at the two end nodes of each element to yield the following two discrete element equations:
Equation (19) can be expressed in a compact matrix indicial form to yield a system of element equations:
Equation (20) provides both the primary variable and its flux at each node of the problem domain and unlike FEM, it does not determine the spatial derivative of the primary variable by an approximation procedure.
3. Results and Discussions
1) Example 1
We test the formulation developed herein by computing the velocity profiles for different values of pressure gradient corresponding to the Falkner-Skan equation. Figure 2 shows that for accelerated flows we obtain velocity profiles without any point of inflection. On the other hand, points of inflection occur for all retarded flows. This is in line with the physics of boundary layer flows. For boundary layer
Figure 2. Velocity profiles for different pressure gradients.
flows across obstacles or blunt bodies, the fluid kinetic energy is dissipated by viscous drag and the fluid is deprived of enough energy to continue along its path. As a result, it halts at some “separation” point after which it reverses and departs from the surface. This is intuitive picture is confirmed by observing the point of inflection in the profile for. Some important examples are considered for the flat plate, and for stagnation point flow. For, the resulting Falkner-Skan equation can be transformed into the differential equation of the rotationally symmetrical flow with stagnation point. For the dependent variable; for (Schlichting  ). This means that we can relate a two-dimensional flow over a wedge to that of a rotationally symmetrical geometry.
2) Example 2
The code developed herein is again validated by solving an example given by Chow  . The energy equation, is addressed by first obtaining expressions from the velocity components from the Blasius equation subject to the three boundary conditions (8). Solutions are obtained by using the Runge-Kutta numerical technique and guessing the starting values of at and then using the iterative technique outlined above. A reasonably large is chosen to approximate infinity. The boundary condition (11b) becomes at. This defines the problem domain for the modified boundary integral domain discretization.
In order to be consistent with Chow  , three cases are considered for the fluid media water, air and mercury with Prandtl numbers 6.750, 0.714 and 0.044 respectively. Table 1 shows the temperature profiles obtained for the different media. These are found to be almost identical with those of Chow  Figure 3 further illustrates the physics of the results. Fluids with higher Prandtl numbers (less thermal conduction) displayed higher temperature profiles. It also confirms the fact that the thermal boundary-layer thickness increases with decreasing Prandtl number due to the increased conductivity. The height is enough to display the profile for water but barely sufficient for air and mercury.
3) Example 3
We solve the Falkner-Skan equation (Equation (6)) together with the energy equation (Equation (10)) to quantify the effect of different pressure gradients on the dimensionless temperature profiles for the three fluid media mentioned in example 2 (Figure 4(a) and Figure 4(b)). Lower temperature profiles are displayed for. This is in consonance with physics. As the fluid approaches separation from a solid surface it looses much of its energy by viscous dissipation, the fluid is left with insufficient energy to proceed further. In an effort to maintain a constant throughput of energy, it reverses and “peels” away from the
Figure 3. Flat plate temperature profiles for different fluid media.
Figure 4. (a) Temperature profiles for separating flow; (b) Temperature profiles for accelerated flow.
surface. Table 2 shows the temperature profiles for which by comparison with Table 1 shows the relative magnitudes of the temperature profiles for the flat plate and for fluid separation. Reasons given above are further confirmed.
Table 1. Flat plate tabulated profiles for different fluid media.
Table 2. Tabulated temperature profile.
4) Example 4
Next we consider the dimensionless energy difference distributions by comparing the total energy content in the free stream and the boundary layer. When a fluid particle originally in the free stream of temperature and speed U becomes isentropically decelerated to zero speed at the surface of the plate, in order to conserve energy, its temperature rises to. This is actually the stagnation temperature of the flow. The total energy per unit mass of the fluid at the wall is the sum of the enthalpy and kinetic energy per unit mass and is the same as that of the free stream of magnitude (Chen  ):
Equation (21) gives the impression that the kinetic energy and by default the free stream velocity is fully recovered. But this is not true for real fluids with non-vanishing kinematic viscosity and thermal conductivity k. In actual fact the kinetic energy is “modified” via the temperature at the wall in order to allow for energy conservation and Equation (21) becomes:
The second part on the right hand side of Equation (22) “recovers” the kinetic energy component of the total energy and is called the recovery factor. Table 1 and Table 2 show that the recovery factor is less than unity for air and mercury whose Prandtl numbers are below unity and is greater than unity for water which possesses a Prandtl number greater than 1. Hence the total fluid energy at the wall may be lower or higher than that in the free stream depending on whether the Prandtl number is greater than or less than unity. The work described herein helps us to gain a handle on the total energy distributions for any of the three cases for different values of pressure gradient. For a semi-infinite flat plate, the total energy is given as:
It can be shown that the dimensionless energy difference between freestream and the wall can be represented as:
Figures 5(a)-(d) confirms that for a fluid with a smaller Prandtl number, where the thermal conductivity is
Figure 5. (a) Dimensionless energy difference; (b) Dimensionless energy difference; (c) Dimensionless energy difference; (d) Dimensionless energy difference.
relatively large, more heat is conducted out of the fluid element than is produced within it by friction, hence the overall energy within its boundary layer becomes lower. This is also combined with the fact, that for a case where the fluid is approaching separation there is a further decrease in the overall energy resulting from a decrease in kinetic energy which further hinders the further advance of the fluid particle. This will in some cases deplete the available boundary layer energy to the extent that it becomes much smaller than that of the free stream and yields a negative value for the dimensionless energy difference.
A robust numerical algorithm involving an element-driven hybridization of the boundary integral theory has been applied to the solution of the energy equation for physical relevant flows of the Falkner-Skan equation. To this author’s knowledge, this is probably one of the few times such an approach is taken for this type of problem. By relying on this hybrid numerical technique, we report highly accurate numerical results without resorting to complex numerical integrations. There are a few issues to emphasize here. First the hybrid procedure permitted Equation (10) to be solved in a 1-D domain. There is no need to seek or device an equivalent analog of the same problem in a 2D domain as would be the case for a classical application of the Boundary Element Method (BEM). Such a transformation will give a false sense of accuracy, because the rigor incurred in order to achieve a BEM domain-reduction is more often than not compensated for by accuracy. The second point is that by dealing directly with such problems without an undue resort to approximations and by discretizing the problem domain we open the door for handling problems involving heterogeneity and nonlinearity (Onyejekwe and Onyejekwe  , Onyejekwe  ). The overall significance of our investigation is that it opens the chapter on the development of highly accurate hybrid boundary integral numerical algorithms for boundary layer flows involving similarity transformations. There is still a paucity of information in this field for boundary layer numerical solutions involving similarity transformations.
 Parand, K., Shahini, M. and Mehdi, D. (2010) Solution of a Laminar Boundary Layer Flow via a Numerical Method. Communications in Nonlinear Science and Numerical Simulation, 15, 360-367.
 Shalini, J. and Choudhary, R. (2015) Effects of MHD on Boundary Layer Flow in Porous Medium Due to Exponentially Shrinking Sheet with Slip. Proccedia Engineering, 127, 1203-1210.
 Bhattacharyya, K., Hayat, T. and Alsaedi, A. (2013) Analytic Solution for Magnetohydrodynamic Boundary Layer Flow of a Casson Fluid over a Stretching/Shrinking Sheet with Wall Mass Transfer. Chinese Physics B, 22, Article ID: 024702.
 Gupta, P.S. and Gupta, A.S. (1977) Heat and Mass Transfer on a Stretching Sheet with Suction or Blowing. The Canadian Journal of Chemical Engineering, 55, 744-746.
 Magyari, E. and Keller, B. (1999) Heat and Mass Transfer in Boundary Layers on an Exponentially Stretching Continuous Surface. Journal of Physics D: Applied Physics, 32, 577-585.
 Parand, K., Shahini, M. and Dehghan, M. (2010) Solution of a Laminar Boundary Layer Flow via a Numerical Method. Communications in Nonlinear Science and Numerical Simulation, 15, 360-367.
 Onyejekwe, O.O. and Onyejekwe, O.N. (2011) Numerical Solutions of One-Phase Classical Stefan Problem Using and Enthalpy Green Element Formulation. Advances in Engineering Software, 42, 743-749.
 Onyejekwe, O.O. (2015) A Direct Implementation of a Modified Boundary Integral Formulation for the Extended Fisher-Kolmogorov Equation. Journal of Applied Mathematics and Physics, 3, 1262-1269.