Eﬀect of temperature dependent viscosity on the thermal instability of two-dimensional stagnation point ﬂow

– This paper presents numerical study of thermal instability of a two-dimensional stagnation point ﬂow when the ﬂuid viscosity is assumed to vary as a linear function of temperature. Similarity transformation was used to reduce the partial diﬀerential boundary layer equations to a non linear system of coupled ordinary diﬀerential equations before solving it numerically using the fourth order Runge-Kutta method with shooting technique. The linear stability of the basic ﬂow to three-dimensional disturbances is then investigated by making use of the normal mode decomposition within the G¨ortler-Hammerlin framework. The equations of linear stability theory create an eigenvalue problem which is solved numerically by means of a pseudo spectral collocation method using Laguerre’s polynomials. The numerical experiment reveals that temperature-dependent viscosity aﬀects signiﬁcantly the onset of thermal instability. It is found that the increase in the temperature-dependent ﬂuid viscosity acts to increase the stability of the basic ﬂow.


Introduction
The two-dimensional stagnation point flow in a mixed convection refers to the flow that results from a twodimensional fluid impinging on a heated surface at right angle and flowing then symmetrically about the stagnation point.This phenomenon is important when the buoyancy forces due to the temperature difference between the surface and the free stream become large.Stagnation flows are found in many applications such as flows over the tips of rockets, aircrafts, submarines and oil ships.In two-dimensional orthogonal stagnation flow, called Hiemenz flow, a similarity transformation is used to reduce the Navier-Stokes equations to a nonlinear ordinary differential equation and leads to a steady solution for the basic flow.A literature survey related to this topic revealed that the stagnation flow problem was started by Hiemenz [1] who developed an exact solution of the Navier-Stokes governing equations for the forced convection case.It is known that the physical properties of fluid may change significantly with temperature.For lubricating fluids, heat generated by the internal friction and the a Corresponding author: mfatsah@hotmail.frcorresponding rise in temperature affects the viscosity of the fluid and so the fluid viscosity can no longer be assumed constant.Therefore, to predict the flow behavior accurately it is necessary to take into account the viscosity variation for incompressible fluids.Gary et al. [2] and Mehta and Sood [3] showed that, when this effect is included the flow characteristics may change substantially compared to the constant viscosity assumption.It is worthy to note that, some other important investigations were conducted by, among others [4][5][6][7][8][9][10][11].In all the aforementioned papers, special emphasis has been given to the effect of variable viscosity on the flow field and heat transfer.As in many other configurations, temperature dependent viscosity affects not only the convective motion of the fluid, but also the stability of the flow.We are not aware of any work for this question in spite of its crucial importance.
In the following, the main concerns were directed to the onset of three-dimensional instability in stagnation point flow.We refer to this flow as the Hiemenz flow.A literature survey related to this topic revealed that Görtler [12] was the first to study the stability problem of steady two-dimensional boundary layer flow of an Article published by EDP Sciences incompressible viscous fluid, and derived the disturbance equations for the Hiemenz flow.These equations were studied by Hammerlin [13] and demonstrated that plane stagnation flow can sustain three-dimensional disturbances.Hiemenz's flow is known to be linearly stable since the work by Wilson and Gladwell [14] who argued that the disturbances must decay exponentially far upstream.Lyell and Huerre [15] confirmed the linear results of Wilson and Gladwell [14].They also presented nonlinear results, obtained with a Galerkin method, suggesting that the flow is unstable for disturbances of large amplitudes.Hall et al. [16] extended the work of Wilson and Gladwell [14] to include both the effect of crossflow in the freestream and the effect of suction or blowing at the wall.Brattkus and Davis [17] considered a wider class of disturbances for the stability of the incompressible Hiemenz flow and concluded that these modes were more stable compared with the self-similar disturbances.
However, to the best of the authors' knowledge, very little is known about the thermal instability of the stagnation point flows over heated horizontal surfaces.In this regard, Chen et al. [18], reveal that thermal excitation generates instabilities when the Rayleigh number exceeds some critical value.It is also found that the critical Rayleigh number and the critical wavenumber are relatively insensitive to the Prandtl number, when defined on the basis of the thermal boundary layer length-scale.These findings remain qualitatively unchanged when obliqueness of the free stream [19] and the relative direction of buoyancy forces [20] are taken into account.Soon thereafter, Amaouche et al. [21] studied the effect of a constant magnetic field on the thermal instability of a two-dimensional stagnation point flow; they found that magnetic fields act to increase its stability.In an experimental investigation carried out by Kitamura and Mitsuishi [22] on mixed convective flows induced around horizontal heated plates keen interests were directed to the occurrence of unstable flows due to buoyancy force.
The present paper will focus on the linear stability of the incompressible stagnation point flow within the Görtler-Hammerlin framework, when the fluid viscosity is assumed to vary as a linear function of temperature.The equations of linear stability theory create an eigenvalue problem which is solved numerically by means of a pseudo spectral collocation method based upon Laguerre's functions.Special attention will be paid to the analysis of the onset conditions of three-dimensional instability and the effects of pertinent parameters such as Prandtl number, Grashof number and viscosity parameter on the linear stability of the basic flow.It's found that the onset conditions of thermal instability are significantly affected by the viscosity parameter and that the increase in the temperature-dependent fluid viscosity acts to increase the stability of the basic flow.

Basic equations
Consider the steady, two-dimensional flow of a viscous and incompressible fluid near the stagnation-point on a heated horizontal flat plate.The geometry of the flow and the coordinate system are indicated in Figure 1.The flat plate is in the (x , z ) plane and the basic flow is twodimensional in the (x , y ) plane, the x axis being the streamwise direction along the plate, and the z axis the spanwise direction.The stagnation streamline coincides with the y axis.The external velocity is prescribed as v e (ax , −ay , 0); where a is constant with a ≥ 0. The flat plate is maintained at a fixed temperature T w higher than the ambient viscous fluid temperature T ∞ .
Here, the physical properties of the fluid are assumed to be constant except the dynamic viscosity which is a linear function of temperature and the density in the buoyancy force term, which is satisfied by the Boussinesq approximation.Under these assumptions, the equations of continuity, momentum, energy and associated boundary conditions are cast in their dimensional form as follow: The boundary conditions for the problem under consideration are In the above equations v is the velocity vector of the fluid, t the time, p the pressure, g the gravitational acceleration and T the temperature.The subscripts w and ∞ stand respectively for the wall and free stream conditions, stars indicate dimensional quantities.ρ, υ, β and α (λ/ρCp) are the density, kinematic viscosity, thermal expansion coefficient, and thermal diffusivity of the fluid, respectively.It is known that for liquid the viscosity decreases as temperature increases.Here the fluid viscosity is assumed to vary as a linear function of temperature as [23]: where μ * is the constant value of the coefficient of viscosity far away from the sheet δ, γ are constants with γ > 0. We choose δ = 1.

Basic flow
The continuity equation is satisfied if we choose a stream function.Equations ( 1)-( 3) can be advantageously rewritten using stream function formalism.Taking the curl of Equation ( 2) allows to eliminate the pressure and the following set of equations is then obtained.
The boundary conditions take the form: The above equations are, for convenience, phrased in terms of dimensionless variables.Distances and time are scaled using the factors = (υ * /a) 1/2 and a −1 , respectively [21].The stream function ψ is referred to υ * ; the scaled temperature is defined by θ The parameters in the problem are Prandtl number P r = υ * /α, viscosity parameter A = γ (T w − T ∞ ) and Grashof number The governing partial differential Equations ( 1)-(3) admit similarity solutions for obtaining the dimensionless stream function f (y) and temperature θ (y).The relative parameters are introduced as: After introducing the similarity transformation, the equations of ( 6) and ( 7) can be transformed into a set of following forms in terms with f (y) and θ (y) that could be expressed as: where the prime denotes a partial differentiation with respect to y.The transformed boundary conditions are given by: The nonlinear differential Equations ( 11) and ( 12) along with the boundary conditions ( 13) and ( 14) are solved numerically using the fourth order Runge-Kutta method with shooting technique.The grid mesh of Δy = 10 −2 is selected to be satisfactory for a convergence criterion of 10 −6 in nearly all cases.

Linear stability analysis
In order to represent a physically pertinent flow, the above basic solution must be stable at least to small three-dimensional disturbances in some subset of parameter space.Therefore we consider infinitesimally small disturbances propagating along the boundary layer, so that the instantaneous quantities ū, v, w, p and θ can be expressed as the sum of the basic-state (u, v, 0, p, θ) and the disturbance-state ũ, ṽ, w, p, θ quantities, so that: ū, v, w, p, θ (x, y, z, t) = (u, v, 0, p, θ) (x, y) + ũ, ṽ, w, p, θ (x, y, z, t) (15) Substituting the above expression into the governing equations of continuity, momentum and energy, subtracting the basic state, and linearizing with respect to the small perturbations gives a set of linearized governing equations of continuity, momentum and energy which best describes the stability characteristics of small perturbations.Substituting the decomposition (15) into Equations ( 1)-( 3), after linearization around the base state we obtain the following system of equations where R 1 and R 2 are the differential operators defined as In this paper, we proceed from the linear governing equations of continuity, momentum and energy equations without making any further simplifications, and recognize that the problem consists of a set of coupled fivedimensional partial-differential equations with variable coefficients.These coefficients, which depend on the basic flow, change strongly in the normal (y) direction and linearly in the chordwise (x) direction, but not in the spanwise (z) direction.As a consequence, the solution is separable in the variables z and t, and retaining the self similarity for the perturbation amplitude, the disturbance quantities of a general travelling mode can be expressed in the form [21] ũ, ṽ, w, p, θ (x, y, z, t) = xû, v, ŵ, p, θ (y) e (ikz+ωt) (21) where k denotes the spanwise wavenumber and ω the temporal growth rate and xû, v, ŵ, p, θ are complex amplitude functions of three-dimensional small disturbances.This particular x-dependence of the two-dimensional disturbance was first proposed and used by Görtler [12] and Hämmerlin [13] in studying the stability of two dimensional stagnation-point flows.These stability results are limited to the special class of self similar disturbances considered above.This so called Görtler-Hammerlin model was extended for a three-dimensional stability analysis by Criminale et al. [24] and Theofilis et al. [25].In the present problem it is customary to focus attention on temporal instability and in general it is assumed that ω is real [26].
Hence, substituting the decomposition (21) into Equations ( 16)-( 20) one obtains The boundary conditions are set as: where The number of unknowns involved in Equations ( 22)-( 26) is reduced by eliminating from the system the pressure and the transverse component of the velocity.This is achieved by using the continuity and transverse momentum equations leading to the following eigenvalue problem A 1 q = ωA 2 q (29) with This system is to be solved numerically in order to determine the eigenvalue ω along with the corresponding eigenfunctions û, v and θ as functions of the spanwise wave number k for the three parameters characterizing the overall fluid system and the base state, Gr, P r and A.
The steady solution is linearly stable if ω > 0 and the disturbance decay with time, otherwise (i.e., ω < 0) it is unstable and the disturbance grows with time.At marginality, corresponding to ω = 0, only for some characteristic values of Gr the problem will allow a non trivial solution for given k, P r and A.

Computational method
In order to approximate the solution of the differential eigenvalue system to be solved and incorporate its exponential damping at infinity, we shall use a pseudo spectral method based on expansion by Laguerre's functions, i.e., Laguerre's polynomials L n (y) multiplied by a decaying exponential [21].The most important feature of the spectral method is exponential convergence which allows for attaining high degree accuracy with a modest number of collocation points.However, the use of Laguerre's polynomials is motivated by the particular distribution of their zeroes, i.e. the first zeroes (used as collocation nodes) are very close together compared to other.This distribution is quite suitable to describe the regions of high gradients like the boundary layer.So, truncating the expansion after a finite number N of terms, an approximation F N of a function F is thought in the form Φ N (y) exp (−y) with Φ N being a polynomial of degree at most N and forced to fulfill linear system at collocation nodes which are selected to be the zeroes of L N (y).Detailed theory of the numerical approach can be found in Reference [27].Taking into account the homogenous boundary condition at y 0 , yields an expression for (û N , vN , θN )(y) of the form: The discretization procedure transforms the original differential eigenvalue problem into an algebraic eigenvalue problem expressed in terms of discretization matrices Q 1 and Q 2 such that with where The N × N square matrices I and D are the identity matrix and the matrix associated to the differential operator D, respectively.Φ N is a vector of expansion coefficients.Hence, as outlined in Reference [21] the problem is posed as a linear generalized eigenvalue problem for ω.To have a nontrivial solution, it is necessary that the matrix Q 1 −ωQ 2 must be singular or, equivalently, its determinant should vanish.This allows computing the ω spectrum for any combination of parameters k, P r, Gr and A. The basic flow is then temporally unstable or stable whether there exists at least one positive eigenvalue or not.Instability occurs only above some level of heating given by the lowest value of Gr, obtained by varying the wavenumber, which cancels det(Q 1 ) at marginality

Results and discussion
The generalized algebraic eigenvalue problem Equation (32) has been solved numerically using a pseudo spectral method based on expansion by Laguerre's functions with implicit QR algorithm.In order to achieve a satisfactory convergence, checks on the influence of the truncation level N were performed in terms of its effects on the basic flow and on the location of the critical conditions for the onset of instability.This showed that the accuracy of the numerical scheme can be improved by increasing the number of collocation nodes.However, higher accuracy is achieved at the expense of a significant increase in computational time.The convergence criteria are based on the relative difference between the current and previous iterative values of critical Grashof number Gr c .When the difference reaches 10 −4 , the solution is assumed to have converged and the iterative process is terminated.
To validate the accuracy of the numerical scheme, we make a comparison of the obtained results corresponding to the neutral stability curve with the available published results of Amaouche and Boukari [20] and Amaouche et al. [21] as shown in Figure 2. As it seen, despite the small deviations observed for small wavenumber k, a good agreement is observed between these two results.These deviations are due to the higher number of collocation nodes N , used in the present study comparatively with above references.Herein after, the numerical computations are performed for several values of viscosity parameter A and Prandtl number Pr.The stability results, overall, depict that the onset conditions of thermal instability are significantly affected by the viscosity parameter.Globally one can observe that viscosity parameter acts to reduce instability, this means that the flow is more stable with variable viscosity.An overview of the stability properties of the basic flow can be seen from the sequence of neutral stability curves displayed in Figure 3 for different values of A and Pr.To generate a neutral curve (ω = 0) we used Newton's method.The iteration process is repeated until |det (Q 1 ) | vanishes within the assumed tolerance (in the present paper |det(Q 1 )| ≤ 10 −6 ).We note that the unstable region lies above the curve while the stable region lies below it.Also each curve shows the minimum value of Grashof number (Gr c ) required to destabilize the mode of wavenumber k.If Gr is below Gr c , the flow will be linearly stable while when Gr exceeds this value, the flow will be unstable.We can also see that basic flow is lin- early unstable when Grashof number exceeds some critical value since the works of Chen et al. [18]; who noted that thermal excitation generates three-dimensional unstable disturbances.It is observed from these figures that the stability of the flow is significantly affected by the viscosity parameter.Specifically, the transition to secondary flow with variable viscosity occurs at a critical Grashof number that may be several times greater than the critical value for classical Hiemenz flow with constant viscosity.We can see that the instability region in the (k, Gr ) plane is always contained in that corresponding to a smaller value of viscosity parameter.This effect is inversed for Pr (i.e.increasing Pr leads to the expansion of the instability region).In addition, it is showed that the stability of the basic flow increases with increasing the viscosity parameter.
Also, Figure 3b shows that the thermal buoyancy variation has no significant influence on the wavenumber of the most unstable modes for smaller values of Pr number (P r < 0.5).In our opinion, the interpretation of these results can be explained physically.Indeed, this can be justified by the fact that, for small values of Pr, thermal disturbances are prone to be dissipated and the most unstable mode stays insensitive to the change on the Gr number.Otherwise, this means that a greater amount of heating is needed to destabilize the basic flow [21].
This observation is confirmed in a global picture provided by Figure 4 showing the critical Grashof number as functions of Prandtl number for different values of viscosity parameter.These figures depict that an increase in Pr leads to reduce critical Grashof number Gr c until a critical value of Prandtl number Pr c corresponding to (Gr c ) min then the effect of increasing Pr leads to increase in Gr c .However, the change in Gr c becomes imperceptible when P r → ∞ Also, we note that the effect of viscosity parameter is imperceptible for small values of Prandtl number Pr.
Results for the least stable and the most unstable branches are plotted in Figure 5 which shows the temporal growth rate ω versus the wavenumber k for different values of viscosity parameter.It is seen that an increase in viscosity parameter leads to the attenuation of the most unstable and the least stable modes.This means that an increase in viscosity parameter acts to retard the transition of the basic flow.This confirms the stabilization effect of viscosity parameter.

Conclusion
In this paper we investigated the stability of steady two-dimensional stagnation point flow when the fluid viscosity is assumed to vary as a linear function of temperature.The eigenvalue problem has been constituted by successfully applying the linear stability theory which is solved numerically by means of a pseudo spectral collocation method using Laguerre's polynomials.It's found that viscosity parameter acts to increase the stability of the basic flow.Furthermore an increase in temperaturedependent fluid viscosity yields (i) more important critical Grashof number for the onset of thermal instability, (ii) the reduction of the region of instability in the (k, Gr ) plane, (iii) the decrease of the temporal growth rate of the most unstable mode and of the least stable one.

Fig. 4 .
Fig. 4. Critical Grashof number Grc versus P r for different values of A.

Fig. 5 .
Fig. 5. Plot of the temporal growth rate versus the wavenumber k.