Incompressible Flow and Heat Transfer over a Plate : A Hybrid Integral Domain-Discretized Numerical Procedure

This work deals with incompressible two-dimensional viscous flow over a semi-infinite plate according to the approximations resulting from Prandtl boundary layer theory. The governing nonlinear coupled partial differential equations describing laminar flow are converted to a self-similar type third order ordinary differential equation known as the Falkner-Skan equation. For the purposes of a numerical solution, the Falkner-Skan equation is converted to a system of first order ordinary differential equations. These are numerically addressed by the conventional shooting and bisection methods coupled with the Runge-Kutta technique. However the accompanying energy equation lends itself to a hybrid numerical finite element-boundary integral application. An appropriate complementary differential equation as well as the Green second identity paves the way for the integral representation of the energy equation. This is followed by a finite element-type discretization of the problem domain. Based on the quality of the results obtained herein, a strong case is made for a hybrid numerical scheme as a useful approach for the numerical resolution of boundary layer flows and species transport. Thanks to the sparsity of the resulting coefficient matrix, the solution profiles not only agree with those of similar problems in literature but also are in consonance with the physics they represent.


Introduction
Ever since its first appearance in literature [1] 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 [1] or its more general counterpart, the Flakner-Skan (F-S) equation [2].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 onedimensional 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 [3]).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 [4] [5].
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 [6], Bhattacharyya et al. [7]).Earlier work involving flow past a stretching plate was recorded by Crane [8].His work was further extended to deal with a permeable surface by Gupta and Gupta [9].Cases involving unsteady shrinking film were considered by Wang [10].This work was later vastly improved by Miklavčič and Wang [11] and Adhikary and Misra [12].Similar work was carried out by Magyari and Keller [13] for self-similar boundary-layer flows induced by permeable stretching walls, as well as for stagnation point flow (Bhattacharyya [14]).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 elementdriven 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.

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 [15]): 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 external flow velocity and the pressure gradient are given by: ( ) where U is the external 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 [2] which are similar but not identical to those of Blasius [1].Non-dimensionalized flow coordinate and stream function and flow velocities yield: A nondimensional stream function ( ) The velocity components become ( ) And A governing equation for f can be found by substituting these non-dimensionalized variables into the momentum equation to yield ( ) where 2 1 Equation ( 6) is the well known Falkner-Skan equation.It comes with the following boundary conditions: and reduces to the Blasius equation for 0 β = .In the case of accelerating and decelerating flows 0, 0 respectively.Most physically realistic flows exist in the region of 0.19884 2 β − ≤ ≤ .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. [16].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 f η ′  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 1 0 at 0 , 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 equ- ation.We non-dimensionalize the energy equation by introducing a dimensionless temperature difference ( ) θ η as: where 1 , T U are the ambient temperature and the free stream velocity respectively, and p c is the specific heat at constant pressure.After substitution into Equation (1c) we obtain: where Pr 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 , i η η are field and source nodes respectively.The first derivative of G with respect to η is , i H x x is the Heaviside function defined as ( ) The resulting integral equation yields: where 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 0 η to L η 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: Equation ( 18) is the interpretation of Equation ( 16) in elemental context where the variable "e" stands for an element and ( ) 1 e η and ( ) 2 e η 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 , and , e e e e θ ϕ θ ϕ .For the numerical implementation of Equation ( 18), we shall evaluate the source node at the two end nodes of each element , e e e e i i η η η η = = 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: M e e e e e ij j ij j i e i j 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.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 η = Ψ (Schlichting [17]).This means that we can relate a two-dimensional flow over a wedge to that of a rotationally symmetrical geometry.

Results and Discussions
2) Example 2 The code developed herein is again validated by solving an example given by Chow [15].The energy equation, is addressed by first obtaining expressions from the velocity components from the Blasius equation .This defines the problem domain for the modified boundary integral domain discretization.In order to be consistent with Chow [15], 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 [15] 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 max 10 η = 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 .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

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 1 T 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 [15]): ( ) 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 ( ) 0 θ 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 ( ) 0 θ 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: 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.

Conclusion
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 [18], Onyejekwe [19]).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.
number and k is the thermal conductivity of the fluid.The transformed boundary conditions are:

1) Example 1
We test the formulation developed herein by computing the velocity profiles ( )f η ′for different values of pressure gradient β corresponding to the Falkner-Skan equation.

Figure 2 Figure 2 .
Figure 2. Velocity profiles for different pressure gradients.

.
Some important examples are considered for the flat plate ( ) into the differential equation of the rotationally symmetrical flow with stagnation point.For the dependent variable ( ) three boundary conditions(8).Solutions are obtained by using the Runge-Kutta numerical technique and guessing the starting values of θ at 0 η = and then using the iterative technique outlined above.A reasonably large max η is chosen to approximate infinity.The boundary condition (11b) becomes 0 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 0.198 β = −

)Figures 5 Figure 5 .
Figures 5(a)-(d) confirms that for a fluid with a smaller Prandtl number, where the thermal conductivity is

Table 2
shows the temperature profiles for 0.198β = −which by comparison with Table1shows the relative magnitudes of the temperature profiles for the flat plate ( )

Table 1 .
Flat plate tabulated profiles for different fluid media ( )