A Combined Analytic and Numerical Procedure for Stagnation-Point Flow

The two dimensional vector form of the Navier-Stokes equation is reduced to a fourth–order equation for the streamfunction. Boundary conditions arise from considerations of the no-slip constraint at the surface as well the interaction of viscous flow with potential-flow at the edge of the boundary layer. By employing a separation of variables technique and introducing certain dimensionless variables, the stream function equation is converted into its dimensionless analog with appropriate boundary conditions. The resulting quasi-linear third-order ordinary differential equation facilitates the numerical computation of the velocity and the pressure terms. This is achieved by solving the nonlinear two-point boundary-value problem with a time-marching method involving a Crank-Nicolson and Newton-linearization schemes until steady-state solution is obtained. The velocity, stream-function and pressure profiles are discussed with reference to various computation parameters and are found to be in good agreement with the physics of the problem. It was also found that there is no penalty in accuracy for a broad range of CFL numbers. However as the CFL number exceeds a certain threshold, the approach to convergence becomes erratic as indicated by the spurious results produced by the solution residuals.


Introduction
It is well known that that the flow at stagnation point, serves as a canonical model of a leading-edge boundary layer. Research activity in this area, closely followed the discovery of boundary-layer concept itself. The measurement of the pressure field on the surface of a cylinder by Prandtl led to the discovery by Hiemenz (1911) of an exact solution of the Navier-Stokes equations which fully describes the stagnation point flow and which up till today bears his name.
Stagnation point flow has long been a special area of research activity for fluid dynamicists because of its certain implications. The heat transfer rates are usually maximum in the stagnation region and are of special interest in areas such as atmospheric reentry of space crafts, in rarefied hypersonic flows and in power output devices. The impulsive interaction of a fluid with a solid surface as well the entrainment of ambient fluid separates it from Blasius flow past a flat plate. Studies involving stagnation point flow have embraced many areas of computational enquiry for example in hydromagnetic flows it was chosen by Na (1979) to illustrate the application of finite difference method (FDM) to a thirdorder boundary value problem. It is also applicable to boundary layers along material handling conveyors, to arterial blood flow problems as well as aerodynamic extrusion of plastic sheets. Later studies have also extended in numerous ways to include various physical considerations. The study of the flow field in stagnation point flow near a stretching sheet in nanofluids is worthy of note. Makinde and Mishra (2015) studied the stagnation point flow of a variable viscosity fluid past a stretching surface with radiative heat. They carried out a parametric study involving the effects of various thermophysical properties on a nanofluid velocity, temperature, concentration, skin friction, Nusselt and Sherwood numbers. A stagnation point flow of MHD chemically reacting nanofluid over a stretching convective surface with slip and radiative heat was carried out by Makinde et al. (2016). They considered the thermophoresis and Brownian motion on hydromagnetic stagnation point flow of a nanofluid with heat and mass transfer over a stretching convective surface. By the use of similarity transformation, the nonlinear governing partial differential was reduced to a set of nonlinear ordinary differential equations which were then solved by the shooting method in conjunction with the Runge-Kutta Fehlberg technique. In another related study Agbaje et al. (2015) adopted a spectral perturbation method (SPM) to study the governing equations of flow and heat transfer of an incompressible electrically conducting fluid near the stagnation point on a stretching sheet. They concluded that the SPM scheme could handle complicated expansion more efficiently than the ordinary perturbation schemes.
Perturbation techniques have also been applied to stagnation point flow for non-Newtonian fluids. Such flows exhibit huge numerical challenges because their constitutive equations generate momentum equations which are multiplied by the fluid viscoelastic parameters. As a result, the order of the governing differential equation exceeds the number of given boundary conditions and the problem converts to a singular perturbation type. The problem was first successfully tackled by Rajeshwari and Rathna (1962) using the Karman-Polhausen method. It was followed by the application of linear perturbation technique by Beard and Walters (1964). Numerical work in this area has been rather slow as seen by subsequent work by Davies (1966), Ng (1981) and Serth (1974). The most recent of these can be traced to Teipel (1986), Ariel (1992). Ariel (1993) introduced a numerical algorithm for computing the stagnation point flow of a second grade fluid with/without suction when the fluid is extracted from the plane at a uniform rate. He obtained asymptotic solutions valid for large values of the suction parameter. A similar attempt involving a perturbation technique was adopted by Massoudi and Ramezan (1992). Their analysis was found to be valid only for small values of the parameter that characterizes non-Newtonian fluids. A later improvement on their work was carried out Garg (1994) who used a pseudo-similarity technique.
Because of the demands of modern technology, research on stagnation point flow is still an ongoing process. Studies have been carried out in such diverse areas as petroleum industries, flow towards solid boundaries in groundwater flows, flow of water towards bridge piers, towards hurricanes, permeable surfaces and geothermal energy extractions, etc. Studies in engineered colloidal particles involving stagnation point flow have also started featuring prominently. Ibrahim and Makinde (2015) examined the effect of slip and convective boundary condition on magnetohydrodynamic (MHD) stagnation point flow and heat transfer due to Casson nanofluid past a stretching sheet and determined that the skin friction coefficient increases with an increase in Casson parameter and decreases with an increase in velocity ratio parameter. In a related work, Ibrahim and Shankar (2013) studied boundary layer flow and heat transfer of a nanofluid over a vertical plate with a convective surface boundary condition. They concluded that the local Nusselt and Sherwood numbers increase with an increase in convective parameter and Nusselt number. We may mention in passing that Casson flow is a prototype non-Newtonian fluid model and relates very much to viscoelastic fluids. We find such models very vital in the study of the flow behavior of industrial products like malt, liquid chocolate, syrup, tomato cream, fruit juices, honey, flow of blood through arteries and so on. Hence most of the earlier work done by Ariel (1993) on stagnation point flow for non-Newtonian fluid as well as those of Massoudi and Ramezan (1992) and Ibrahim and Makinde (2015) should be considered as benchmark work in this area of research.
In the work reported herein, we deployed the vorticity equation to convert the two-dimensional, second order vector form of the Navier-Stokes equation to a fourth-order differential equation for the stream-function. The accompanying boundary conditions were accounted for by invoking the no-slip constraint at the surface as well as the constraint imposed by the inviscid flow regime at the far field. Adopting the separation of variables technique and introducing some dimensionless variables, the fourth order governing differential equation together with the boundary conditions were converted into a desired dimensionless thirdorder differential equation. This equation originates from the Navier-Stokes equation and can be regarded as the viscousflow replacement of the Laplace's equation for the streamfunction. Having established the relationship between the stream-function and the horizontal velocity, a time marching numerical procedure was introduced to solve the governing equation to steady state. Once we compute the velocity field, we can determine the pressure field by taking the divergence of the momentum equation to subsequently arrive at the Poisson equation for a two-dimensional incompressible flow.

Model Formulation
Consider the laminar two-dimensional flow an incompressible fluid impinging on a plate situated at y=0, ( Figure 1).
where u, v are components of the velocity vector. The continuity equation for the two-dimensional potential flow is given by: Stagnation point flow is not vorticity free and is proportional to the Laplacian of the stream-function.
The vorticity is transported from the wall by the dual processes of convection and diffusion and the appropriate vorticity transport equation for a two-dimensional flow is given as: where ν is the kinematic viscosity.
Substituting equation (3) into equation (4) yields the equation for the stream-function.
where 4 ∇ is the biharmonic operator. We hasten to comment on a few interesting features concerning equation (5). In addition to the continuity equation being unconditionally satisfied, this single equation combined with appropriate boundary conditions adequately describes a two-dimensional fluid flow and avoids some of the numerical complications inherent in the Navier-Stokes equation. It has only the kinematic viscosity ν as a parameter. We also note that equation for creeping flow is satisfied when the left side is equal to zero. From equation (5), the required fourth order partial differential equation for the stream-function is given by.
Appropriate boundary conditions are assigned to equation (6); two of them are due to the no-slip requirement at the surface and the second is assigned as a result of the interaction between the viscid and inviscid regions of the flow field far away from the plate.
We note in passing that since the boundary layer thickness is not known a priori, the boundary condition at the edge of the boundary layer is specified at a far enough distance from the plate that is at infinity so to speak. Typically one would solve the governing equation or the stream function ψ with associated boundary conditions, instead of solving the Navier-Stokes and continuity equations directly for the velocity and pressure fields. Once the stream function is determined, It is possible to compute the velocity and pressure components. The price to be paid for this procedure is that the stream-function satisfies a fourth-order partial differential equation in contrast to the to the Navier-Stokes equation which is a second order partial differential equation. To enhance the numerical computation of equation (6) we first of all consider certain aspects of its analytic features. Since the potential velocity v depends on y and the other component u depends on x, the following form of stream function is suggested as solution.
Substituting equation (8) in to equation (6) reveals the following details about the components of the stream function ( ) In addition, ( ) y Φ satisfies the following equation.
with accompanying boundary conditions given as: We can now recast the stream-function equations in terms of the following dimensionless variables in line with previous work by Ariel ( The link between stream-function and velocity is established by: Equations (10a) and (10b) can be rewritten as: Equation (12) has many attractive features, it obviates a rather restrictive approach of scaling down equation (10a) into systems of ordinary differential equations and corresponding boundary conditions. It admits and facilitates the application of many domain and boundary based techniques to the numerical solution of stagnation point flows. In the work reported herein, we have elected to apply a finite difference time-marching method to proceed to a steady state numerical solution and then determine the stream function by a trapezoidal rule resolution of the following Once we have determined the stream-function and the velocity fields, then it is time to handle the pressure. In order to achieve the pressure Poisson equation formulation, we take the divergence of the momentum equation and use the divergence free condition. We initiate this procedure by considering the x and y components of the 2-D Navier-Stokes equation, differentiate the first with respect to x and carry out the same procedure for the second with respect to y, add both and rearrange to finally yield: Equation (13) is the Poisson equation for pressure which ensures that continuity is satisfied. It can be seen that both the pressure and velocity field are now coupled in a continuous domain. We still need to relate appropriate boundary conditions to equation (13) otherwise the computation will be fraught with errors in the pressure field and in the incompressibility conditions and the proposed scheme can not produce faithful results all the way to the boundary. For this paper, we employ a different approach in a similar spirit to the one used by Wilcox (ibid) and Papanastasiou (ibid). We observe that for stagnation point flow, Bernoulli's equation can sometimes lead to some very interesting conclusions regarding flow along a streamline. And for such cases, the stagnation pressure is the sum of static and dynamic pressures. Because of the dynamic link between the pressure and velocity fields we expect both to approach the inviscid solution in the limit as y → ∞ ( Figure  1). This idea can be deployed in an easily computable form by noting that.
where stag p is the stagnation pressure, and the second term of the right of equation (14) is the dynamic pressure. The following solution is suggested.
Substituting equations (15a) and (10a) into equation (13) yields the following ( ) 2 X x x = and: ( ) y as y dy dy dy Equation (15b) is an easily computable analog of equation (13). And the boundary conditions comply with the requirement that the pressure distribution agree with the inviscid value far from the location of the stagnation point. Following the procedure adopted in this study equation (15b) is put in a dimensionless form by introducing some dimensionless variables The resulting second order ODE is represented as: We carry out a further simplification of equation (16) Finally an explicit equation for the pressure field is obtained by integrating equation (17) twice and applying the boundary conditions for equation (16) to yield.
The overall physics of the formulation is now fully justified. The stream-function is directly related to the velocity field whose changes result in acceleration of the fluid particles. At the same time, local accelerations are controlled by the vorticity field while the pressure field adjusts accordingly to establish a balance of the net force due to pressure and viscosity.

Problem Discretization
We recast the two-point boundary value problem (equation 12 a) in the following manner.

( )
The boundary conditions are the same as specified in equation (12b). For the initial condition we assume a quadratic variation of ( ) u η near the surface. The discretization scheme we shall adopt here is the Crank-Nicolson-Newton-Richtmeyer scheme. To this end equation The following stencils are employed to handle the right side of equation (19b) By the same token: Similarly the second component of the second derivative reads: The general expression for a nonlinear term Equation (19) is put together to form a tridiagonal matrix in terms of the new dependent variable ( ) and then solved with the Thomas algorithm. The above scheme can be applied straightforwardly to quantify the simulation residuals: After some numerical experimentation, we situated the upper finite-difference grid at 5 η = with the hope that it satisfies the condition at infinity. In order to enhance stability, we have adopted the time step limitation and the effective CFL number set by Anderson (1990):

Results and Discusion
The algorithm developed herein is tested on some typical examples for stagnation-point flows. The interconnectedness between our three variables of interest namely velocity, stream-function and pressure has very well been established. Table 1 shows the numerical results for the horizontal velocity, pressure and stream-function with respect to the dimensionless independent variable η . The outcome in each case not only agrees with the specified boundary conditions but also displays the unique relationships between them. We go a step further to display this by plotting the three variables with respect to the independent variable η . Figure 2 confirms that the horizontal velocity is zero at the surface ( ) 0 η = according to the noslip boundary constraint and asymptotically approaches the free-stream value far away from the stagnation.   (11).
Using the boundary condition ( ) 0 0 ϕ = , the rest of the profile can be evaluated. In addition, Figure 2 shows that unlike the horizontal velocity profile, the profile of the stream function does not approach the free-stream asymptotically though they both start at zero at the surface. The shape of the streamfunction profile in Figure 3 can be ascertained by looking intuitively at the area under the curve of the velocity profile that results from the relationship of the two dependent In addition the boundary conditions are satisfied at both ends of the profile. Once the stream-function is determined, calculation of the pressure field near the stagnation point is facilitated. This is the effective stagnation pressure in the inviscid flow above the stagnation point. It has both static and dynamic components. As with the velocity, it approaches the far-field inviscid solution as the vertical dimensionless coordinate η approaches infinity. Figure 4 is in complete agreement with equation (16) which determines the trajectory of the pressure field and is different from the actual stagnation pressure due to contribution from its dynamic component.    η ≤ ≤ the velocity field displays an asymptotic profile. This accounts for the descending part of the residual profile. Some 'wiggles' are observed in the region around 2. 5 3.1 η ≤ ≤ which corresponds roughly with the edge of the boundary layer profile.

Conclusions
The 2-dimensional stagnation point flow of an incompressible flow has been studied with a Crank-Nicolson-Newton-Richtmeyer time scheme. Stable results were obtained with 51 grid points and a CFL number of 25. Plots produced herein in complete agreement with boundary layer theory. A noticeable aspect of this study is that it not only facilitates the computation of the pressure field which is hardly mentioned in many papers on stagnation point flow, the time matching procedure obviates the need for the introduction of systems of ODE for the eventual solution of the boundary-layer type governing equations. As a consequence it paves the way for the application of boundary and domain based numerical techniques to such problems. [8] Garg, V. K., Heat transfer due to stagnation point flow on a non-Newtonian fluid, Acta Mech. Vol. 104 159-171 (1994).