**Computation of Non-Isothermal Reversed Stagnation-Point Flow over a Flat Plate**

Vai Kuong Sin and Chon Kit Chio *University of Macau Taipa, Macao SAR, China*

#### **1. Introduction**

158 Computational Simulations and Applications

Zaichik, L. I. & Alipchenkov, V. M. (2005). Statistical Models for Predicting Particle

Thermal Calculation of Power Generators (Standard Method). Moscow, Energija, 1973, 295

*of Heat and Fluid Flow*, Vol.26, pp. 416-430.

pages [in Russian].

Dispersion and Preferential Concentration in Turbulent Flows. *International Journal* 

The full Navier-Stokes equations are difficult or impossible to obtain an exact solution in almost every real situation because of the analytic difficulties associated with the nonlinearity due to convective acceleration. The existence of exact solutions are fundamental not only in their own right as solutions of particular flows, but also are agreeable in accuracy checks for numerical solutions.

In some simplified cases, such as two-dimensional stagnation point flows, by introducing coordinate variable transformation, the number of independent variables is reduced by one or more. The governing equations can be simplified to the non-linear ordinary differential equations and are analytic solvable. The classic problems of two-dimensional stagnation-point flows can be analyzed exactly by Hiemenz Hiemenz (1911), one of Prandtl's first students. These are exact solutions for flow directed perpendicular to an infinite flat plate. Howarth Howarth (1951) and Davey Davey (1961) extended the two-dimensional and axisymmetric flows to three dimensions, which is based on boundary layer approximation in the direction normal to the plane.

The similarity solutions for the temperature field were studied by Eckert Eckert (1942). Case corresponding a step change in wall temperature or in wall heat flux in laminar steady flows at a stagnation point has been also investigated by several authors (see Chao et al. Chao & Jeng (1965), Sano Sano (1981) and Gorla Gorla (1988)). Further, Lok et al. Lok et al. (2006) investigated the mixed convection near non-orthogonal stagnation point flow on a vertical plate with uniform surface heat flux, where the results published are very good with present value of *θ*(0) for the constant wall temperature boundary condition. On the contrary, reversed stagnation-point flow over an infinite flat wall does not have analytic solution in two-dimensional steady state case.

The aim of this study is to investigate the unsteady viscous reversed non-isothermal stagnation-point flow, which is started impulsively in motion with a constant velocity away from near the stagnation point. A similarity solution of full Navier-Stokes equations and energy equation are solved by applying numerical method. Studies of the reversed stagnation-point flow have been considered during the last few years, as this flow can be applied in different important applications that occur in oil recovery industry, as shown in Fig. (1).

*ψ* = −

*∂u <sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>u</sup>*

*∂v <sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>u</sup>*

concentrate on the transient behavior in other section.

temperature distribution in the reversed stagnation-point flow.

From the definiton of stream function, we have

equation

variable *η*.

**2.2 Energy equation**

with the boundary conditions

*∂u ∂x* + *v ∂u <sup>∂</sup><sup>y</sup>* <sup>=</sup> <sup>−</sup><sup>1</sup> *ρ ∂p ∂x* + *ν*( *∂*2*u <sup>∂</sup>x*<sup>2</sup> <sup>+</sup>

*∂v ∂x* + *v ∂v <sup>∂</sup><sup>y</sup>* <sup>=</sup> <sup>−</sup><sup>1</sup> *ρ ∂p ∂y* + *ν*( *∂*2*v <sup>∂</sup>x*<sup>2</sup> <sup>+</sup>

*η* =

where *u* and *v* are the velocity components along *x* and *y* axes, and *ρ* is the density.

*<sup>u</sup>* <sup>=</sup> *∂ψ*

*<sup>v</sup>* <sup>=</sup> <sup>−</sup>*∂ψ*

Substituting *u* and *v* into the governing equations results a simplified partial differential

Equation (6) is the similarity equation of the full Navier-Stokes equations at two-dimension reversed stagnation point. The coordinates *x* and *y* are vanished, leaving only a dimensionless

When the flow is in steady state such that *fητ* = 0, it can be proved that the differential equation does not have solution under the boundary conditions Davey (1961). Thus we

In this section, we shall focus on the non-isothermal flow which is at a temperature *T* different from that of the wall *Tw*. By solving the energy equation, we are able to determine the

*A*

Computation of Non-Isothermal Reversed Stagnation-Point Flow over a Flat Plate 161

where *η* is the non-dimensional distance from wall and *τ* is the non-dimensional time. Noting that the stream function automatically satisfies equation of continuity 1 . The Navier-Stokes equations White (2003) governing the unsteady flow with constant physical properties are

<sup>√</sup>*Aνx f*(*η*, *<sup>τ</sup>*) (3a)

*τ* = *At* (3c)

*∂*2*u*

*∂*2*v*

*<sup>∂</sup><sup>y</sup>* <sup>=</sup> <sup>−</sup>*Axf<sup>η</sup>* (5a)

*<sup>∂</sup><sup>x</sup>* <sup>=</sup> <sup>√</sup>*A<sup>ν</sup> <sup>f</sup>* (5b)

*<sup>f</sup>ητ* <sup>−</sup> (*fη*)<sup>2</sup> <sup>+</sup> *f fηη* <sup>−</sup> *<sup>f</sup>ηηη* <sup>+</sup> <sup>1</sup> <sup>=</sup> 0, (6)

*f*(0, *τ*) = *fη*(0, *τ*) = 0 (7a) *fη*(∞, *τ*) = 1. (7b)

*<sup>ν</sup> <sup>y</sup>* (3b)

*<sup>∂</sup>y*<sup>2</sup> ) (4a)

*<sup>∂</sup>y*<sup>2</sup> ) (4b)

Fig. 1. Oil recovery industry

#### **2. Governing equation**

#### **2.1 Momentum equation**

The viscous fluid flows in a rectangular Cartesian coordinates (*x*, *y*, *z*), Fig. 2, which illustrates the motion of external flow directly moves perpendicular out of an infinite flat plane wall. The origin is the so-called stagnation point and *z* is the normal to the plane.

Fig. 2. Cooridnate system of reversed stagnation-point flow

By conservation of mass principle with constant physical properties , the equation of continuity is:

$$\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 \tag{1}$$

We consider the two-dimensional reversed stagnation-point flow in unsteady state and the flow is bounded by an infinite plane *y* = 0, the fluid remains at rest when time *t* < 0. At *t* = 0, it starts impulsively in motion which is determined by the stream function

$$
\psi = -axy \tag{2}
$$

For a viscid fluid the stream function, since the flow motion is determined by only two factors, the kinematic viscosity *ν* and *α*, we consider the following modified stream function

$$
\psi = -\sqrt{A\nu}\mathfrak{x}f(\eta, \mathfrak{r})\tag{3a}
$$

$$
\eta = \sqrt{\frac{A}{\nu}} y
$$

$$
\boldsymbol{\pi} = \boldsymbol{A}\boldsymbol{t} \tag{3c}
$$

where *η* is the non-dimensional distance from wall and *τ* is the non-dimensional time. Noting that the stream function automatically satisfies equation of continuity 1 . The Navier-Stokes equations White (2003) governing the unsteady flow with constant physical properties are

$$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial \mathbf{x}} + v \frac{\partial u}{\partial y} = -\frac{1}{\rho} \frac{\partial p}{\partial \mathbf{x}} + \nu (\frac{\partial^2 u}{\partial \mathbf{x}^2} + \frac{\partial^2 u}{\partial y^2}) \tag{4a}$$

$$\frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} = -\frac{1}{\rho} \frac{\partial p}{\partial y} + \nu (\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}) \tag{4b}$$

where *u* and *v* are the velocity components along *x* and *y* axes, and *ρ* is the density. From the definiton of stream function, we have

$$
\mu = \frac{\partial \psi}{\partial y} = -A x f\_{\eta} \tag{5a}
$$

$$v = -\frac{\partial \psi}{\partial x} = \sqrt{A\nu}f\tag{5b}$$

Substituting *u* and *v* into the governing equations results a simplified partial differential equation

$$(f\_{\eta\tau} - (f\_{\eta})^2 + ff\_{\eta\eta} - f\_{\eta\eta\eta} + 1 = 0,\tag{6}$$

with the boundary conditions

2 Will-be-set-by-IN-TECH

The viscous fluid flows in a rectangular Cartesian coordinates (*x*, *y*, *z*), Fig. 2, which illustrates the motion of external flow directly moves perpendicular out of an infinite flat plane wall. The

*y*

✻

✕ ✕ ❑ ❑

✻

By conservation of mass principle with constant physical properties , the equation of

We consider the two-dimensional reversed stagnation-point flow in unsteady state and the flow is bounded by an infinite plane *y* = 0, the fluid remains at rest when time *t* < 0. At *t* = 0,

For a viscid fluid the stream function, since the flow motion is determined by only two factors,

the kinematic viscosity *ν* and *α*, we consider the following modified stream function

*u* = 0 *O*

*∂u ∂x* + *∂v*

it starts impulsively in motion which is determined by the stream function

Fig. 2. Cooridnate system of reversed stagnation-point flow

✲

*u* = −*Ax*

*<sup>∂</sup><sup>y</sup>* <sup>=</sup> <sup>0</sup> (1)

*ψ* = −*αxy* (2)

*x*

origin is the so-called stagnation point and *z* is the normal to the plane.

Fig. 1. Oil recovery industry

**2. Governing equation**

**2.1 Momentum equation**

continuity is:

$$f(0, \tau) = f\_{\eta}(0, \tau) = 0 \tag{7a}$$

$$f\_{\parallel}(\infty,\tau) = 1.\tag{7b}$$

Equation (6) is the similarity equation of the full Navier-Stokes equations at two-dimension reversed stagnation point. The coordinates *x* and *y* are vanished, leaving only a dimensionless variable *η*.

When the flow is in steady state such that *fητ* = 0, it can be proved that the differential equation does not have solution under the boundary conditions Davey (1961). Thus we concentrate on the transient behavior in other section.

#### **2.2 Energy equation**

In this section, we shall focus on the non-isothermal flow which is at a temperature *T* different from that of the wall *Tw*. By solving the energy equation, we are able to determine the temperature distribution in the reversed stagnation-point flow.

boundary. On the contrary, the viscous forces were neglected away from the wall. The convection terms dominate the motion of external flow in considering the inviscid equation in the fluid. According to their solution, the general features of the predicted streamline are

Computation of Non-Isothermal Reversed Stagnation-Point Flow over a Flat Plate 163

*<sup>f</sup>ητ* <sup>−</sup> (*fη*)<sup>2</sup> <sup>+</sup> *f fηη* <sup>+</sup> <sup>1</sup> <sup>=</sup> 0. (15)

*f* = −*η*, *f* � = −1 (18)

)<sup>2</sup> <sup>−</sup> <sup>1</sup> <sup>=</sup> <sup>0</sup> (19)

(0) = 0 (20a)

(∞) = −1 (20b)

*θ*(*η*, *τ*) = Θ(*γ*) (21)

*<sup>τ</sup>F*(*γ*) (16)

<sup>−</sup>*cγ*) (17)

We therefore consider the similarity of the inviscid equation

and the further integration provides an exact solution

Proudman and Johnson obtained a similarity solution of (15) is in the form

*f*(*η*, *τ*) = *e*

*c* (1 − *e*

where c is a constant of integration; the improved numerical evaluations of Robins and Howarth Robins & Howarth (1972) estimate the value of *c* to be approximately 3.51. This solution describes the flow in the outer region. When *<sup>τ</sup>* <sup>→</sup> <sup>∞</sup> and *<sup>η</sup>*/*e<sup>τ</sup>* is relatively small, the

*F* ∼ −*γ* = −*ηe*

*f* ��� − *f f* �� + (*f* �

*f*(0) = *f* �

This is exactly the classic stagnation-point problem Hiemenz (1911)) by changing the sign in *f* . It is a third-order nonlinear ordinary differential equation and does not have an analytic solution, and thus it is necessary to solve it numerically and the result is show in Fig. 4 .

Wen *τ* → ∞, we have an exact solution of the momentum equation in outer region, and thus, we still apply the same procedure to solve the temperature profiles in outer region. Consider

> *γ* = *ηe* −*τ*

*f* �

−*τ*

*<sup>F</sup>*(*γ*) = *<sup>γ</sup>* <sup>−</sup> <sup>2</sup>

sketched in Fig. (3).

solution (17) yields

Substituting in equation (6) becomes

with the boundary conditions

**3.2 Temperature distribution**

the following transformation:

where

and

For constant-property fluid such as results, the transient energy equation Burmeister (1993) is given as follow

$$
\rho c\_p \left( \frac{\partial T}{\partial t} + \mu \frac{\partial T}{\partial x} + v \frac{\partial T}{\partial y} \right) = k \left( \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} \right) + \mu \Phi \tag{8}
$$

where *k* is the thermal conductivity and *cp* is the heat capacity. Note that *u* and *v* are the velocity components. These give

$$
\mu = -A \mathbf{x} f\_{\eta} \tag{9a}
$$

$$v = \sqrt{Av}f\tag{9b}$$

and Φ is defined as

$$\Phi = 2\left[\left(\frac{\partial u}{\partial \mathbf{x}}\right)^2 + \left(\frac{\partial v}{\partial y}\right)^2\right] + \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial \mathbf{x}}\right)^2 - \frac{2}{3}\left(\frac{\partial u}{\partial \mathbf{x}} + \frac{\partial v}{\partial y}\right)^2\tag{10}$$

and is called the viscous dissipation since it represents the irreversible conservation of mechanical forms of energy to a thermal form.

To transform equation (8) into a nondimensional form, a dimensionless temperature *θ* is defined as

$$\theta = \frac{T - T\_w}{T\_{\infty} - T\_w}.\tag{11}$$

where *T*<sup>∞</sup> is the ambient temperature. Noting that both *Tw* and *T*<sup>∞</sup> are constant,

$$T(0,t) = T\_{\overline{w}\_{\prime}} \ T(\infty,t) = T\_{\infty \prime} \tag{12}$$

the temperature distribution *T* can be considered as a function of *y* and *t* only. Under the assumption that the viscous dissipation is negligible compared to conduction at the wall, the energy equation may be written as

$$
\theta\_{\eta\eta} - \frac{\rho c\_p \nu}{k} f \theta\_\eta = \frac{\rho c\_p \nu}{k} \theta\_\tau \tag{13}
$$

subject to the boundary conditions

$$
\theta(0,\tau) = 0 \qquad \theta(\infty,\tau) = 1 \tag{14}
$$

Equation (13) is a second-order partial differential equation with variable coefficients *f*(*η*, *τ*) and the Prandtl number *Pr* = *ρcpν*/*k* is assumed to be constant. Consider the fluid of which *Pr* = 1, the thermal boundary layer and the velocity boundary layer collapse, and thus, substituting *θ* = *f* � , equation (6) and (13) represent the same equation.

#### **3. Asymptotic solution**

#### **3.1 Velocity distribution**

When *τ* is relatively small, Proudman and Johnson Proudman & Johnson (1962) first considered the early stages of the diffusion of the initial vortex sheet at *y* = 0. They suggested that, when the flow is near the plane region, the viscous forces are dominant, and the viscous term in the governing Navier-Stokes equations is important only near the boundary. On the contrary, the viscous forces were neglected away from the wall. The convection terms dominate the motion of external flow in considering the inviscid equation in the fluid. According to their solution, the general features of the predicted streamline are sketched in Fig. (3).

We therefore consider the similarity of the inviscid equation

$$\left(f\_{\eta\tau} - \left(f\_{\eta}\right)^2 + f f\_{\eta\eta} + 1 = 0.\tag{15}$$

Proudman and Johnson obtained a similarity solution of (15) is in the form

$$f(\eta, \tau) = e^{\tau} F(\gamma) \tag{16}$$

and the further integration provides an exact solution

$$F(\gamma) = \gamma - \frac{2}{c}(1 - e^{-c\gamma})\tag{17}$$

where c is a constant of integration; the improved numerical evaluations of Robins and Howarth Robins & Howarth (1972) estimate the value of *c* to be approximately 3.51. This solution describes the flow in the outer region. When *<sup>τ</sup>* <sup>→</sup> <sup>∞</sup> and *<sup>η</sup>*/*e<sup>τ</sup>* is relatively small, the solution (17) yields

−*τ*

*F* ∼ −*γ* = −*ηe*

and

4 Will-be-set-by-IN-TECH

For constant-property fluid such as results, the transient energy equation Burmeister (1993) is

and is called the viscous dissipation since it represents the irreversible conservation of

To transform equation (8) into a nondimensional form, a dimensionless temperature *θ* is

*<sup>θ</sup>* <sup>=</sup> *<sup>T</sup>* <sup>−</sup> *Tw T*<sup>∞</sup> − *Tw*

Under the assumption that the viscous dissipation is negligible compared to conduction at the

Equation (13) is a second-order partial differential equation with variable coefficients *f*(*η*, *τ*) and the Prandtl number *Pr* = *ρcpν*/*k* is assumed to be constant. Consider the fluid of which *Pr* = 1, the thermal boundary layer and the velocity boundary layer collapse, and thus,

When *τ* is relatively small, Proudman and Johnson Proudman & Johnson (1962) first considered the early stages of the diffusion of the initial vortex sheet at *y* = 0. They suggested that, when the flow is near the plane region, the viscous forces are dominant, and the viscous term in the governing Navier-Stokes equations is important only near the

, equation (6) and (13) represent the same equation.

*<sup>k</sup> <sup>f</sup> θη* <sup>=</sup> *<sup>ρ</sup>cp<sup>ν</sup>*

where *T*<sup>∞</sup> is the ambient temperature. Noting that both *Tw* and *T*<sup>∞</sup> are constant,

the temperature distribution *T* can be considered as a function of *y* and *t* only.

*θηη* <sup>−</sup> *<sup>ρ</sup>cp<sup>ν</sup>*

 *∂*2*T <sup>∂</sup>x*<sup>2</sup> <sup>+</sup> *∂*2*T ∂y*<sup>2</sup> 

*u* = −*Axf<sup>η</sup>* (9a) *<sup>v</sup>* <sup>=</sup> <sup>√</sup>*A<sup>ν</sup> <sup>f</sup>* (9b)

> *∂u ∂x* + *∂v ∂y* <sup>2</sup>

*T*(0, *t*) = *Tw*, *T*(∞, *t*) = *T*∞, (12)

*θ*(0, *τ*) = 0 *θ*(∞, *τ*) = 1 (14)

. (11)

*<sup>k</sup> θτ* (13)

+ *μ*Φ (8)

(10)

given as follow

and Φ is defined as

defined as

*ρcp*

Φ = 2

 *<sup>∂</sup><sup>u</sup> ∂x* <sup>2</sup> + *∂v ∂y* <sup>2</sup> + *∂u ∂y* + *∂v ∂x* <sup>2</sup> − 2 3

mechanical forms of energy to a thermal form.

wall, the energy equation may be written as

subject to the boundary conditions

substituting *θ* = *f* �

**3. Asymptotic solution**

**3.1 Velocity distribution**

 *∂T <sup>∂</sup><sup>t</sup>* <sup>+</sup> *<sup>u</sup>*

where *k* is the thermal conductivity and *cp* is the heat capacity. Note that *u* and *v* are the velocity components. These give

*∂T ∂x* + *v ∂T ∂y* = *k*

$$f = -\eta, \quad f' = -1 \tag{18}$$

Substituting in equation (6) becomes

$$f'''' - ff'' + (f')^2 - 1 = 0\tag{19}$$

with the boundary conditions

$$f(0) = f'(0) = 0\tag{20a}$$

$$f'(\infty) = -1\tag{20b}$$

This is exactly the classic stagnation-point problem Hiemenz (1911)) by changing the sign in *f* . It is a third-order nonlinear ordinary differential equation and does not have an analytic solution, and thus it is necessary to solve it numerically and the result is show in Fig. 4 .

#### **3.2 Temperature distribution**

Wen *τ* → ∞, we have an exact solution of the momentum equation in outer region, and thus, we still apply the same procedure to solve the temperature profiles in outer region. Consider the following transformation:

$$\theta(\eta,\tau) = \Theta(\gamma) \tag{21}$$

where

$$\gamma = \eta e^{-\tau}$$

if we consider a finite value of *τ*, equation (13) becomes an ordinary differential equation

Computation of Non-Isothermal Reversed Stagnation-Point Flow over a Flat Plate 165

*<sup>c</sup>* (<sup>1</sup> <sup>−</sup> *<sup>e</sup>*

*c* = 3.51. The following pages (Figs. (5) to (6)) show the asymptotic solution of temperature

According to the previous work, the governing equations in reversed stagnation-point flow

The above equations (6) and (13) subject to the boundary conditions are are nonlinear third-order partial differential equations. They do not admit similarity solution and numerical

We shall, however, use here a numerical method. It is an implicit finite-difference method with second-order accuracy. The partial differential equations can be expressed as approximate expressions, so that it is easy to program the solution of large numbers of coupled equation.

<sup>−</sup>*cγ*)Θ� = 0 (22)

Θ(0) = 0, Θ(∞) = 1 (23)

*<sup>f</sup>ητ* <sup>−</sup> (*fη*)<sup>2</sup> <sup>+</sup> *f fηη* <sup>−</sup> *<sup>f</sup>ηηη* <sup>+</sup> <sup>1</sup> <sup>=</sup> 0. (6)

*<sup>f</sup>ητ* <sup>=</sup> *<sup>f</sup>ηηη* + (*fη*)<sup>2</sup> <sup>−</sup> <sup>1</sup> <sup>+</sup> *f fηη* (24a)

*Pr θηη* <sup>−</sup> *<sup>f</sup> θη* (24b)

*h* = 1 − *f<sup>η</sup>* (25a) *g* = *θ* (25b)

(1 − *h*) *dη* (26a)

(1 − *h*) *dη* (26b)

*θηη* − *Pr f θη* = *Pr θτ* (13)

2*Pr e*2*<sup>τ</sup>*

Θ�� +

or perturbation methods are needed to solve the equation.

and introducing the new dependent variables

The equations can be rewritten as

We start with rewriting the partial differential equations in the form:

*θτ* <sup>=</sup> <sup>1</sup>

*<sup>h</sup><sup>τ</sup>* <sup>=</sup> *<sup>h</sup>ηη* <sup>+</sup> <sup>2</sup>*<sup>h</sup>* <sup>−</sup> *<sup>h</sup>*<sup>2</sup> <sup>+</sup> *<sup>h</sup><sup>η</sup>*

*Pr <sup>g</sup>ηη* <sup>−</sup> *<sup>g</sup><sup>η</sup>*

*<sup>g</sup><sup>τ</sup>* <sup>=</sup> <sup>1</sup>

subject to the boundary

**4.1 Governing equations**

distributions with *Pr* in outer region.

**4. Finite-difference formulations**

where

are

Fig. 3. Streamlines of reversed stagnation-point flow

Fig. 4. Numerical solutions of stagnation-point flow

if we consider a finite value of *τ*, equation (13) becomes an ordinary differential equation

$$2\Theta'' + \frac{2\operatorname{Pr}\,e^{2\tau}}{c}(1 - e^{-c\gamma})\Theta' = 0\tag{22}$$

subject to the boundary

$$
\Theta(0) = 0, \ \Theta(\infty) = 1 \tag{23}
$$

where

6 Will-be-set-by-IN-TECH

*y*

✻

✕ ✕ ❑ ❑

✲ ✛

Fig. 3. Streamlines of reversed stagnation-point flow

Fig. 4. Numerical solutions of stagnation-point flow

✲

✛ ✲

*x*

*c* = 3.51.

The following pages (Figs. (5) to (6)) show the asymptotic solution of temperature distributions with *Pr* in outer region.

#### **4. Finite-difference formulations**

#### **4.1 Governing equations**

According to the previous work, the governing equations in reversed stagnation-point flow are

$$(f\_{\eta\tau} - (f\_{\eta})^2 + ff\_{\eta\eta} - f\_{\eta\eta\eta} + 1 = 0. \tag{6}$$

$$
\theta\_{\eta\eta} - \Pr\, f\theta\_{\eta} = \Pr\, \theta\_{\tau} \tag{13}
$$

The above equations (6) and (13) subject to the boundary conditions are are nonlinear third-order partial differential equations. They do not admit similarity solution and numerical or perturbation methods are needed to solve the equation.

We shall, however, use here a numerical method. It is an implicit finite-difference method with second-order accuracy. The partial differential equations can be expressed as approximate expressions, so that it is easy to program the solution of large numbers of coupled equation.

We start with rewriting the partial differential equations in the form:

$$f\_{\eta\tau} = f\_{\eta\eta\eta} + \left(f\_{\eta}\right)^2 - 1 + f f\_{\eta\eta} \tag{24a}$$

$$
\theta\_{\tau} = \frac{1}{Pr} \theta\_{\eta \eta} - f \theta\_{\eta} \tag{24b}
$$

and introducing the new dependent variables

$$h = 1 - f\_{\eta} \tag{25a}$$

$$\mathbf{g} = \theta \,\tag{25b}$$

The equations can be rewritten as

$$h\_{\pi} = h\_{\eta\eta} + 2h - h^2 + h\_{\eta} \int (1 - h) \, d\eta \tag{26a}$$

$$g\_{\tau} = \frac{1}{Pr} g\_{\eta\eta} - g\_{\eta} \int (1 - h) \, d\eta \tag{26b}$$

(a) *Pr* = 3

Computation of Non-Isothermal Reversed Stagnation-Point Flow over a Flat Plate 167

(b) *Pr* = 10

Fig. 6. Asymptotic temperature solution Θ for various value of *γ*

Fig. 5. Asymptotic temperature solution Θ for various value of *γ*

8 Will-be-set-by-IN-TECH

(a) *Pr* = 0.7

(b) *Pr* = 1

Fig. 5. Asymptotic temperature solution Θ for various value of *γ*

Fig. 6. Asymptotic temperature solution Θ for various value of *γ*

The finite-difference form of the ODE is written for the midpoint (*τn*, *ηj*), the discretized

Computation of Non-Isothermal Reversed Stagnation-Point Flow over a Flat Plate 169

*<sup>i</sup>* <sup>+</sup> *<sup>h</sup>n*+<sup>1</sup> *i*−1 (Δ*η*)<sup>2</sup> <sup>+</sup> <sup>2</sup>*h<sup>n</sup>*

> *i*Δ*η* 0

*<sup>i</sup>* <sup>+</sup> *<sup>g</sup>n*+<sup>1</sup> *i*−1

> *i*Δ*η* 0

 2*h<sup>n</sup> <sup>i</sup>* <sup>−</sup> (*h<sup>n</sup>*

*Pr <sup>g</sup>n*+<sup>1</sup> *<sup>i</sup>*−<sup>1</sup> <sup>=</sup> *<sup>g</sup><sup>n</sup>* *<sup>i</sup>* <sup>−</sup> (*h<sup>n</sup> i* )2

*<sup>i</sup>* )<sup>2</sup> <sup>−</sup> *<sup>h</sup><sup>n</sup>*

*<sup>i</sup>* − Δ*τ*

(1 − *h*) *dη* (30a)

(1 − *h*) *dη* (30b)

*i* ∑ 0

> *i* ∑ 0

(<sup>1</sup> <sup>−</sup> *<sup>h</sup><sup>n</sup> i* ) (31a)

(<sup>1</sup> <sup>−</sup> *<sup>h</sup><sup>n</sup>*

*<sup>i</sup>* )(31b)

(33a)

(33b)

*<sup>i</sup>*+<sup>1</sup> <sup>−</sup> *<sup>h</sup><sup>n</sup> i*−1 2

*gn <sup>i</sup>*+<sup>1</sup> <sup>−</sup> *<sup>g</sup><sup>n</sup> i*−1 2

*h<sup>τ</sup>* = *hηη* (32a)

*Pr <sup>g</sup>ηη* (32b)

*exp*(−*ξ*2) *<sup>d</sup><sup>ξ</sup>* (34)

*<sup>i</sup>*+<sup>1</sup> <sup>−</sup> <sup>2</sup>*hn*+<sup>1</sup>

*<sup>i</sup>*+<sup>1</sup> <sup>−</sup> <sup>2</sup>*gn*+<sup>1</sup>

*<sup>i</sup>* + Δ*τ*

*<sup>i</sup>* <sup>−</sup> *<sup>β</sup>*

The initial condition is the solution of the following linear partial differential equations

*<sup>g</sup><sup>τ</sup>* <sup>=</sup> <sup>1</sup>

*<sup>h</sup>* <sup>=</sup> <sup>1</sup> <sup>−</sup> *er f <sup>η</sup>*

*<sup>g</sup>* <sup>=</sup> *er f <sup>η</sup>* 2

<sup>√</sup>*<sup>π</sup>*

 *z* 0

*er f*(*z*) = <sup>2</sup>

2 √*τ* 

<sup>√</sup>*τ*/*Pr*

*Pr*(Δ*η*)<sup>2</sup>

*<sup>i</sup>*+<sup>1</sup> <sup>−</sup> *<sup>h</sup><sup>n</sup> i*−1 2Δ*η*

*<sup>i</sup>*+<sup>1</sup> <sup>−</sup> *<sup>g</sup><sup>n</sup> i*−1 2Δ*η*

− *hn*

− *gn*

*<sup>i</sup>*−<sup>1</sup> <sup>=</sup> *<sup>h</sup><sup>n</sup>*

2*β Pr*)*gn*+<sup>1</sup>

This procedure yields the following linear tridiagonal system:

*<sup>i</sup>* <sup>−</sup> *<sup>β</sup>hn*+<sup>1</sup>

*<sup>i</sup>*+<sup>1</sup> + (1 +

equation takes the form.

<sup>−</sup>*βhn*+<sup>1</sup>

**4.2 Initial flow**

where *β* = Δ*τ*/(Δ*η*)2.

an the required solutions are

where the error function *er f*(*z*) is defined as

*<sup>i</sup>*+<sup>1</sup> + (<sup>1</sup> <sup>+</sup> <sup>2</sup>*β*)*hn*+<sup>1</sup>

− *β Pr <sup>g</sup>n*+<sup>1</sup>

*hn*+<sup>1</sup> *<sup>i</sup>* <sup>−</sup> *<sup>h</sup><sup>n</sup> i* <sup>Δ</sup>*<sup>τ</sup>* <sup>=</sup> *<sup>h</sup>n*+<sup>1</sup>

*gn*+<sup>1</sup> *<sup>i</sup>* <sup>−</sup> *<sup>g</sup><sup>n</sup> i* <sup>Δ</sup>*<sup>τ</sup>* <sup>=</sup> *<sup>g</sup>n*+<sup>1</sup>

We now consider the net rectangle in the *τ* − *η* plane shown in Fig. (7) and the net points defined as below:

$$\begin{array}{ll} \eta^{0} = 0, & \eta\_{j} = \eta\_{j-1} + \Delta \eta, & j = 1,2,... \text{J}, \ \eta\_{J} = \eta\_{\infty} \\\\ \tau^{0} = 0, & \tau^{n} = \tau^{n-1} + + \Delta \tau, & n = 1,2,... \text{J}, \end{array}$$

Here *n* and *j* are just the sequence of numbers that indicate the coordinate location, not tensor indices or exponents. The partial differential equations are easily discretized by central

Fig. 7. Net rectangle for finite-difference method

difference representations with second-order accuracy, for example the finite-difference forms for any points are

$$\hbar\_{\eta} = \frac{\hbar\_{i+1}^n - \hbar\_{i-1}^n}{2\Delta\eta} \tag{27}$$

and

$$h\_{\eta\eta} = \frac{h\_{i+1}^n - 2h\_i^n + h\_{i-1}^n}{\Delta\eta} \tag{28}$$

When *i* = 0, since the value of *h<sup>n</sup> <sup>i</sup>*−<sup>1</sup> is not logical, the derivative is replaced by the forward difference with second-order accuracy

$$h\_{\eta} = \frac{-h\_{i+2}^n + 4h\_{i+1}^n - 3h\_i^n}{2\Delta\eta} \tag{29}$$

The finite-difference form of the ODE is written for the midpoint (*τn*, *ηj*), the discretized equation takes the form.

$$\frac{h\_i^{n+1} - h\_i^n}{\Delta \tau} = \frac{h\_{i+1}^{n+1} - 2h\_i^{n+1} + h\_{i-1}^{n+1}}{(\Delta \eta)^2} + 2h\_i^n - (h\_i^n)^2$$

$$-\frac{h\_{i+1}^n - h\_{i-1}^n}{2\Delta \eta} \int\_0^{i\Delta \eta} (1 - h) \, d\eta \tag{30a}$$

$$\frac{\mathcal{S}\_i^{n+1} - \mathcal{S}\_i^n}{\Delta \tau} = \frac{\mathcal{S}\_{i+1}^{n+1} - 2\mathcal{S}\_i^{n+1} + \mathcal{S}\_{i-1}^{n+1}}{Pr(\Delta \eta)^2}$$

$$-\frac{\mathcal{S}\_{i+1}^{\eta} - \mathcal{S}\_{i-1}^{\eta}}{2\Delta \eta} \int\_0^{i\Delta \eta} (1 - h) \, d\eta \tag{30b}$$

This procedure yields the following linear tridiagonal system:

$$-\beta h\_{i+1}^{n+1} + (1+2\beta)h\_i^{n+1} - \beta h\_{i-1}^{n+1} = h\_i^n + \Delta \tau \left[ 2h\_i^n - (h\_i^n)^2 - \frac{h\_{i+1}^n - h\_{i-1}^n}{2} \sum\_{0}^i (1 - h\_i^n) \right] \text{ (31a)}$$

$$-\frac{\beta}{Pr} g\_{i+1}^{n+1} + (1 + \frac{2\beta}{Pr}) g\_i^{n+1} - \frac{\beta}{Pr} g\_{i-1}^{n+1} = g\_i^n - \Delta \tau \frac{g\_{i+1}^n - g\_{i-1}^n}{2} \sum\_{0}^i (1 - h\_i^n) \text{ (31b)}$$

where *β* = Δ*τ*/(Δ*η*)2.

#### **4.2 Initial flow**

10 Will-be-set-by-IN-TECH

We now consider the net rectangle in the *τ* − *η* plane shown in Fig. (7) and the net points

*<sup>η</sup>*<sup>0</sup> <sup>=</sup> 0, *<sup>η</sup><sup>j</sup>* <sup>=</sup> *<sup>η</sup>j*−<sup>1</sup> <sup>+</sup> <sup>Δ</sup>*η*, *<sup>j</sup>* <sup>=</sup> 1, 2, ...*J*, *<sup>η</sup><sup>J</sup>* <sup>=</sup> *<sup>η</sup>*<sup>∞</sup>

Here *n* and *j* are just the sequence of numbers that indicate the coordinate location, not tensor indices or exponents. The partial differential equations are easily discretized by central

difference representations with second-order accuracy, for example the finite-difference forms

*<sup>i</sup>*+<sup>1</sup> <sup>−</sup> <sup>2</sup>*h<sup>n</sup>*

*<sup>i</sup>*+<sup>2</sup> <sup>+</sup> <sup>4</sup>*h<sup>n</sup>*

*<sup>i</sup>*+<sup>1</sup> <sup>−</sup> *<sup>h</sup><sup>n</sup> i*−1

> *<sup>i</sup>* <sup>+</sup> *<sup>h</sup><sup>n</sup> i*−1

*<sup>i</sup>*+<sup>1</sup> <sup>−</sup> <sup>3</sup>*h<sup>n</sup> i*

<sup>2</sup>Δ*<sup>η</sup>* (27)

<sup>Δ</sup>*<sup>η</sup>* (28)

<sup>2</sup>Δ*<sup>η</sup>* (29)

*<sup>i</sup>*−<sup>1</sup> is not logical, the derivative is replaced by the forward

*<sup>h</sup><sup>η</sup>* <sup>=</sup> *<sup>h</sup><sup>n</sup>*

*<sup>h</sup>ηη* <sup>=</sup> *<sup>h</sup><sup>n</sup>*

*<sup>h</sup><sup>η</sup>* <sup>=</sup> <sup>−</sup>*h<sup>n</sup>*

*τ*<sup>0</sup> = 0, *τ<sup>n</sup>* = *τn*−<sup>1</sup> + +Δ*τ*, *n* = 1, 2, ...*J*,

defined as below:

Fig. 7. Net rectangle for finite-difference method

for any points are

When *i* = 0, since the value of *h<sup>n</sup>*

difference with second-order accuracy

and

The initial condition is the solution of the following linear partial differential equations

$$h\_{\pi} = h\_{\eta \eta} \tag{32a}$$

$$\mathbf{g}\_{\pi} = \frac{1}{Pr} \mathbf{g}\_{\eta \eta} \tag{32b}$$

an the required solutions are

$$h = 1 - \operatorname{erf}\left(\frac{\eta}{2\sqrt{\tau}}\right) \tag{33a}$$

$$\mathbf{g} = \text{erf}\left(\frac{\eta}{2\sqrt{\pi/Pr}}\right) \tag{33b}$$

where the error function *er f*(*z*) is defined as

$$\operatorname{erf}(z) = \frac{2}{\sqrt{\pi}} \int\_0^z \exp(-\tilde{\xi}^2) \, d\tilde{\xi} \tag{34}$$

When *τ* → 0, the boundary conditions are convenient to write in the form

$$h\_0^n = \mathbf{g}\_0^n = \mathbf{0},$$

$$h\_i^0 = \operatorname{erf}\left(\frac{\eta\_i}{2\sqrt{\tau}}\right).$$

$$\mathbf{g}\_i^0 = \operatorname{erf}\left(\frac{\eta\_i}{2\sqrt{\tau/Pr}}\right) \tag{35}$$

Equations (31) are defined as being implicit, as more than one unknown appears in the left hand side. They are unconditionally stable, however, set of linear algebraic equations is required to be solved by the tridiagonal matrix algorithm (TDMA), also know as the Thomas algorithm, which is a simplified form of Gaussian elimination that is applied to evaluate tridiagonal systems of equations.

The procedure is straightforward, except for the algebra. The resulting algorithm of the finite-difference method is written in MATLAB, a numerical computing environment allowing matrix manipulations and plotting of functions and data.

#### **5. Numerical result**

The following figure (Figs. (8)) shows the velocity distribution *f* . From the solution 17, we have

$$
\log(1 - f\_{\eta}) = -c\eta + \log 2\tag{36}
$$

Fig. 8. Numerical solution of Eq. (17) against (a) *τ*, (b) *ηe*−*<sup>τ</sup>*

Computation of Non-Isothermal Reversed Stagnation-Point Flow over a Flat Plate 171

and if the similarity solution holds, then the graph of *log*(<sup>1</sup> <sup>−</sup> *<sup>f</sup>η*) against *<sup>η</sup>e*−*<sup>τ</sup>* should provide a straight line of gradient −*c*, except for small values of *η*, In Figs. (8), the values of *log*(1 − *fη*) are plotted against *ηe*−*τ*−3.5 at different value of *τ*. The value of *c* is 3.51 which is consistent with the previous discussion.

Next, we show the numerical solution of temperature distributions with *Pr* in Figs. (9) through (10)) . It is noticed that the dimensionless wall temperature gradient Θ� (0) raises with increase of Prandtl number, but the thermal boundary layer thickness decrease with increase of Prandtl number. The thermal boundary layer thickness is the distance from the body at which the temperature is 99% of the temperature obtained from an inviscid solution. The decrease of thickness can be explained by the definition of Prandtl number that Prandtl number is inversely proportional to the molecular thermal diffusivity *α*. If the Prandtl number is greater than 1, the thermal boundary layer is thinner than the velocity boundary layer. If the Prandtl number is less than 1, which is the case for air at standard conditions, the thermal boundary layer is thicker than the velocity boundary layer.

In comparison to the asymptotic solution, we note that increase in non-dimensional time *τ* leads to an increase in temperature profiles in both cases. Near the wall region where *γ* is small, the dimensionless wall temperature gradient of the numerical solution is lower than that of the asymptotic solution. It is because the asymptotic solution is only valid for the outer region. At our level of discretization, however, we are only able to resolve in small time range. 12 Will-be-set-by-IN-TECH

<sup>0</sup> = 0,

 *η<sup>i</sup>* 2 √*τ* ,

 *η<sup>i</sup>* 2 <sup>√</sup>*τ*/*Pr*

Equations (31) are defined as being implicit, as more than one unknown appears in the left hand side. They are unconditionally stable, however, set of linear algebraic equations is required to be solved by the tridiagonal matrix algorithm (TDMA), also know as the Thomas algorithm, which is a simplified form of Gaussian elimination that is applied to evaluate

The procedure is straightforward, except for the algebra. The resulting algorithm of the finite-difference method is written in MATLAB, a numerical computing environment allowing

The following figure (Figs. (8)) shows the velocity distribution *f* . From the solution 17, we

and if the similarity solution holds, then the graph of *log*(<sup>1</sup> <sup>−</sup> *<sup>f</sup>η*) against *<sup>η</sup>e*−*<sup>τ</sup>* should provide a straight line of gradient −*c*, except for small values of *η*, In Figs. (8), the values of *log*(1 − *fη*) are plotted against *ηe*−*τ*−3.5 at different value of *τ*. The value of *c* is 3.51 which is consistent

Next, we show the numerical solution of temperature distributions with *Pr* in Figs. (9)

with increase of Prandtl number, but the thermal boundary layer thickness decrease with increase of Prandtl number. The thermal boundary layer thickness is the distance from the body at which the temperature is 99% of the temperature obtained from an inviscid solution. The decrease of thickness can be explained by the definition of Prandtl number that Prandtl number is inversely proportional to the molecular thermal diffusivity *α*. If the Prandtl number is greater than 1, the thermal boundary layer is thinner than the velocity boundary layer. If the Prandtl number is less than 1, which is the case for air at standard conditions, the thermal

In comparison to the asymptotic solution, we note that increase in non-dimensional time *τ* leads to an increase in temperature profiles in both cases. Near the wall region where *γ* is small, the dimensionless wall temperature gradient of the numerical solution is lower than that of the asymptotic solution. It is because the asymptotic solution is only valid for the outer region. At our level of discretization, however, we are only able to resolve in small time range.

through (10)) . It is noticed that the dimensionless wall temperature gradient Θ�

*log*(1 − *fη*) = −*cη* + *log* 2 (36)

(35)

(0) raises

When *τ* → 0, the boundary conditions are convenient to write in the form *hn* <sup>0</sup> = *<sup>g</sup><sup>n</sup>*

> *h*0 *<sup>i</sup>* = *er f*

> *g*0 *<sup>i</sup>* = *er f*

tridiagonal systems of equations.

**5. Numerical result**

with the previous discussion.

have

matrix manipulations and plotting of functions and data.

boundary layer is thicker than the velocity boundary layer.

Fig. 8. Numerical solution of Eq. (17) against (a) *τ*, (b) *ηe*−*<sup>τ</sup>*

(a) *Pr* = 3

Computation of Non-Isothermal Reversed Stagnation-Point Flow over a Flat Plate 173

(b) *Pr* = 10

Fig. 10. Numerical temperature solution Θ for various value of *γ*

Fig. 9. Numerical temperature solution Θ for various value of *γ*

14 Will-be-set-by-IN-TECH

(a) *Pr* = 0.7

(b) *Pr* = 1

Fig. 9. Numerical temperature solution Θ for various value of *γ*

Fig. 10. Numerical temperature solution Θ for various value of *γ*

**Part 2** 

**Applications Involving**

**Computational Fluid Dynamics** 
