Journal of Applied Mathematics
Volume 2008 (2008), Article ID 835380, 16 pages
doi:10.1155/2008/835380
Research Article

On the Asymptotic Approach to Thermosolutal Convection in Heated Slow Reactive Boundary Layer Flows

1Department of Mathematics and Applied Mathematics, University of Venda, Private Bag x5050, Thohoyandou 0950, South Africa
2School of Mathematical Sciences, University of KwaZulu Natal, Private Bag X01, Scottsville, Pietermaritzburg 3209, South Africa
3Mathematics Department, University of Swaziland, Private Bag 4, Kwaluseni, Swaziland

Received 22 April 2008; Revised 14 July 2008; Accepted 12 August 2008

Academic Editor: Jacek Rokicki

Copyright © 2008 Stanford Shateyi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

The study sought to investigate thermosolutal convection and stability of two dimensional disturbances imposed on a heated boundary layer flow over a semi-infinite horizontal plate composed of a chemical species using a self-consistent asymptotic method. The chemical species reacts as it diffuses into the nearby fluid causing density stratification and inducing a buoyancy force. The existence of significant temperature gradients near the plate surface results in additional buoyancy and decrease in viscosity. We derive the linear neutral results by analyzing asymptotically the multideck structure of the perturbed flow in the limit of large Reynolds numbers. The study shows that for small Damkohler numbers, increasing buoyancy has a destabilizing effect on the upper branch Tollmien-Schlichting (TS) instability waves. Similarly, increasing the Damkohler numbers (which corresponds to increasing the reaction rate) has a destabilizing effect on the TS wave modes. However, for small Damkohler numbers, negative buoyancy stabilizes the boundary layer flow.

1. Introduction

Convection in which the buoyancy forces are due to both temperature and chemical concentration gradients are referred to as thermosolutal or double diffusive convection. Ostrach [1] pointed out that different modes of such convection exist depending on how the temperature and concentration gradients are oriented relative to one another. Some natural convection flows in the atmosphere and micrometeorological phenomena are often thermosolutal. The heating of the earth by the sun causes atmospheric thermal convection which is usually modified by the presence of moisture evaporated from the ground.

In lakes and oceans, thermosolutal convection is caused by stable vertical concentration distribution with heating from the side or from the top. The stability theory has been used to explain the occurrence of layered structures observed in oceans as explained by Turner [2].

Studies on natural convection flows caused by the simultaneous diffusion of thermal energy and chemical species were carried out by Gebhart and Pera [3]. They considered small species concentration levels and showed that the species Boussinesq approximation led to similarity solutions similar in form to those for single buoyancy mechanism flows.

The effect of heat transfer on the upper-branch stability of Tollmien-Schlichting instabilities (TSIs) in accelerating boundary layer over a rigid surface in incompressible flows was investigated by Mureithi et al. [4]. The study indicated that buoyancy has destabilizing effect on rigid bodies. Their analysis also showed that in the presence of strong buoyancy forces, the five-zone asymptotic structure alters but for moderate buoyancy, the five-zone structure of Smith and Bodonyi [5] remains with very few alterations.

Shateyi et al. [6] considered the effect of fluid buoyancy and chemical reaction between the chemical species and the fluid on the linear stability of two dimensional disturbances wave modes. They extended the theory of boundary layer flows over horizontal surfaces to include a chemical species and the effect of the Damkohler number. Results showed that when the wave number and speed number are varied against the scaled Damkohler number, the effect of increasing buoyancy was destabilizing in agreement with assertions by Motsa et al. [7]. It was shown that increasing reaction kinematics (thus, reducing the density and viscosity) has destabilizing effects on the TS waves.

Pons and Le Quéré [8] showed that Boussinesq equations do not exactly represent buoyancy—induced natural convection. The study observed that thermodynamic consistency is retrieved when both the work due pressure forces and the heat generated by viscous friction are accounted for in the heat equation. The study recommended that any theoretical study of buoyancy-induced natural convection should be done with the thermodynamic Boussinesq model and not with the usual Boussinesq approximations. However, in spite of this weakness, the usual Boussinesq equations are still very useful (see, e.g., Azizi et al. [9]). To that end in this study, we will use the Boussinesq approximations.

The present work presents an asymptotic analysis of the flow induced by buoyancy effects due to both temperature and chemical concentration gradients near the flat surface. This paper is a direct extension of the earlier work in Shateyi et al. [6] to include the effects of temperature differences within the boundary-layer. The flow has uniform surface conditions with the buoyancy effects primarily away from the surface. Our analysis is limited to processes which occur at low concentration gradients. We give an asymptotic investigation of the interactions between the reaction kinematics and the fluid hydrodynamics with the Damkohler number (the ratio of the hydrodynamic time scale to the reaction time scale) as the parameter of primary interest. In the limit of large Damkohler numbers, the reaction kinematics proceed at a much faster rate compared to the fluid hydrodynamics. If the Damkohler number is close to zero, the chemical reactions are slow compared to the motion of the fluid. In this case, a nonreactive fluid can be assumed. The greatest interaction between reaction kinematics and fluid dynamics occurs when the Damkohler is of 𝑂 ( 1 ) .

In the study, we focus attention on the case when the Damkohler number is small (that is, less than unity). The case for large Damkohler numbers was investigated by Shateyi et al. [6]. The aim is to determine the influence of small Damkohler numbers on the stability characteristics of the upper-branch TS instability waves. The major difference between the current work and that of Mureithi et al. [4] arises from the introduction of reaction kinematics. The absence of wall compliance and the presence of a chemical species makes the current work different from that of Motsa et al. [7]. The presence of wall heating or cooling and small Damkohler numbers makes the current work different from Shateyi et al. [6]. The current approach provides asymptotic solutions to the linearised Navier-Stokes rather the numerical solution of the Orr-Sommerfeld equation.

2. Mathematical Formulation

We consider a two-dimensional, incompressible fluid flow over a heated horizontal plate which is composed of a chemical species maintained at a fixed concentration. The chemical species diffuses into the nearby fluid inducing a buoyancy force. A change in the temperature of the fluid near the plate surface due to exogenous heating effects also results in additional buoyancy.

The governing nondimensional unsteady Navier-Stokes equations for an incompressible fluid under a Boussinesq-type approximation are as follows: 𝜕 𝑢 + 𝜕 𝑥 𝜕 𝑣 𝜕 𝑦 = 0 , ( 2 . 1 ) 𝜕 𝑢 𝜕 𝑡 + 𝑢 𝜕 𝑢 𝜕 𝑥 + 𝑣 𝜕 𝑢 𝜕 𝑦 = 𝜕 𝑝 + 1 𝜕 𝑥 ( 𝜕 R e 2 𝑢 𝜕 𝑥 2 + 𝜕 2 𝑢 𝜕 𝑦 2 ) , ( 2 . 2 ) 𝜕 𝑣 𝜕 𝑡 + 𝑢 𝜕 𝑣 𝜕 𝑥 + 𝑣 𝜕 𝑣 𝜕 𝑦 = 𝜕 𝑝 + 1 𝜕 𝑦 ( 𝜕 R e 2 𝑣 𝜕 𝑥 2 + 𝜕 2 𝑣 𝜕 𝑦 2 ) + 𝐺 𝑡 𝑇 + 𝐺 𝑐 𝐶 , ( 2 . 3 ) 𝜕 𝐶 𝜕 𝑡 + 𝑢 𝜕 𝐶 𝜕 𝑥 + 𝑣 𝜕 𝐶 = 1 𝜕 𝑦 ( 𝜕 R e S c 2 𝐶 𝜕 𝑥 2 + 𝜕 2 𝐶 𝜕 𝑦 2 ) D a 𝐶 , ( 2 . 4 ) 𝜕 𝑇 𝜕 𝑡 + 𝑢 𝜕 𝑇 𝜕 𝑥 + 𝑣 𝜕 𝑇 = 1 𝜕 𝑦 ( 𝜕 P r R e 2 𝑇 𝜕 𝑥 2 + 𝜕 2 𝑇 𝜕 𝑦 2 ) . ( 2 . 5 ) However, (2.1)–(2.5) have been nondimensionalised such that the space coordinates are given by ( 𝑥 , 𝑦 ) = 𝐿 ( 𝑥 , 𝑦 ) , the velocity components are ( 𝑢 , 𝑣 ) = 𝑈 ( 𝑢 , 𝑣 ) , the pressure is 𝑝 = 𝜌 𝑈 𝑝 , and the time is 𝑡 = ( 𝐿 / 𝑈 ) 𝑡 , where 𝐿 is the characteristic length scale (e.g., the distance measured from the leading edge of the plate). The species chemical concentration and the temperature are, respectively, 𝐶 𝐶 = 𝐶 𝐶 𝑤 𝐶 𝑇 , 𝑇 = 𝑇 𝑇 𝑤 𝑇 , ( 2 . 6 ) where the asterisks in these relations denote dimensional quantities and the subscripts “ 𝑤 " and “ " refer to the conditions at the wall and the free stream values, respectively. For purposes of this study, the important parameters are the Damkohler number D a that is defined as the ratio of the flow time scale to the chemical time scale and the buoyancy terms 𝐺 𝑐 = 𝐺 𝑟 𝑐 / R e 2 and 𝐺 𝑡 = 𝐺 𝑟 𝑡 / R e 2 , where 𝐺 𝑟 𝑐 = 𝛽 𝑐 𝑔 𝐿 3 ( 𝐶 𝑤 𝐶 ) / 𝜈 2 and 𝐺 𝑟 𝑡 = 𝛽 𝑡 𝑔 𝐿 3 ( 𝑇 𝑤 𝑇 ) are, respectively, the Grashof numbers for mass and thermal diffusion, 𝛽 𝑐 is the coefficient of expansion with respect to mass transfer, and 𝛽 𝑡 is the volumetric coefficient of thermal expansion

The other parameters and variables in (2.1)–(2.5) are the Prandtl number P r , the flow Reynolds number R e , and the Schmidt number S c ( = 𝜇 / 𝐷 𝜌 ) , where 𝐷 is the binary diffusion coefficient, 𝜇 is the coefficient of dynamic viscosity, 𝜌 is the density of the fluid, 𝜈 is the coefficient of kinematic viscosity, and 𝑔 is the acceleration due to gravity.

The momentum boundary conditions are the no-slip conditions: 𝑢 = 𝑣 = 0 a t 𝑦 = 0 . ( 2 . 7 ) We assume that the horizontal plate is maintained at a prescribed temperature 𝑇 𝐵 𝑤 and uniform concentration 𝐶 𝐵 𝑤 . In the far field, we assume that the fluid temperature, species chemical concentration, and the fluid velocity approach their free stream values.

The basic boundary-layer flow is given by 𝑈 ( 𝑢 , 𝑣 , 𝑝 , 𝑇 , 𝐶 ) = 𝐵 , R e 1 / 2 𝑉 𝐵 , 𝑝 𝐵 , 𝑇 𝐵 , 𝐶 𝐵 ( 𝑥 , 𝑌 ) , ( 2 . 8 ) where 𝑌 = R e 1 / 2 𝑦 is the boundary coordinate. For general pressure gradient boundary layers, the following properties hold: 𝑈 𝐵 𝜆 1 𝑌 + 𝜆 2 𝑌 2 𝑇 + a s 𝑌 0 , 𝐵 𝛾 1 + 𝛾 2 𝑈 𝑌 + a s 𝑌 0 , 𝐵 𝜇 1 + 𝜇 2 𝑈 𝑌 + a s 𝑌 0 , 𝐵 1 , 𝑇 𝐵 0 , 𝐶 𝐵 0 a s 𝑌 . ( 2 . 9 ) The coefficients 𝜆 1 = 𝑈 𝐵 𝑦 | 𝑦 = 0 > 0 and 𝜆 2 = 𝑈 𝐵 𝑦 𝑦 | 𝑦 = 0 < 0 are, respectively, the skin friction and curvature of the basic flow profile. The coefficients 𝛾 1 and 𝛾 2 are the heat transfer coefficients and 𝜇 1 and 𝜇 2 are the concentration transfer coefficients.

Disturbances to the basic flow of amplitude factor 𝛿 are introduced and these spread through the boundary layer. The Reynolds number is large while the amplitude of the spatially growing disturbances is assumed to be sufficiently small for linear theory to hold. The TSI grows by extracting energy from the mean flow to the disturbance within the boundary layer (see, e.g., Carpenter and Gajjar [10]).

The asymptotic structure of the boundary layer is now well known and described (e.g., Smith and Bodonyi [5] and Mureithi et al. [4]). It consists of the five regions (see Figure 1) where 𝑅 1 is the main part of the boundary layer with thickness 𝑂 ( R e ( 1 / 2 ) ) , 𝑅 2 is a thinner inviscid adjustment region of thickness 𝑂 ( R e ( 7 / 1 2 ) ) containing the viscous critical layer 𝑅 3 , the wall layer 𝑅 4 of thickness 𝑂 ( R e ( 2 / 3 ) ) , and the outer potential flow 𝑅 5 of thickness 𝑂 ( R e ( 5 / 1 2 ) ) . Consistent with earlier work on linear analysis, the viscous critical layer is ignored in the present theory except in so far as it produces a phase shift in the boundary layer pressure and velocity.

835380.fig.001
Figure 1: Schematic sketch of flow structure showing the multilayered nature of the boundary layer and the relative positioning of the five regions.

The disturbances are taken to be in the form of a modulated wave train periodic in 𝑋 where, for the upper branch of the neutral stability curve, the scaled streamwise and temporal variables are 𝑥 = 𝜖 5 𝑋 , 𝑡 = 𝜖 4 𝜏 , and 𝜖 = R e ( 1 / 1 2 ) is a small parameter. The neutral wave number 𝛼 and phase speed 𝑐 of the disturbances are 𝑂 ( 𝜖 5 ) and 𝑂 ( 𝜖 4 ) , respectively. We thus expand 𝛼 and 𝑐 as 𝛼 = 𝜖 5 𝛼 0 + , 𝑐 = 𝜖 4 𝑐 0 + , ( 2 . 1 0 ) where 𝛼 0 and 𝑐 0 are the scaled real wave number and real phase speed of the travelling wave disturbances. The derivatives 𝜕 / 𝜕 𝑥 and 𝜕 / 𝜕 𝑡 are then replaced by 𝜖 5 𝛼 0 𝜕 / 𝜕 𝑋 and 𝜖 4 𝜔 0 𝜕 / 𝜕 𝜏 , respectively. The other important scalings (see also Shateyi et al. [6]) are 𝐺 𝑐 / 𝐺 𝑡 = 𝑂 ( 𝜖 5 ) and D a 𝑂 ( 𝜖 4 ) when the TS eigenrelation is significantly altered for the first time by the effects of fluid buoyancy.

3. Stability Analysis

Information on the disturbance expansions relevant to the upper-branch stability of boundary layers is now well documented in the literature (see, e.g., Gajjar and Smith [11], Motsa et al. [7], and the references therein). Only the details necessary to obtain the linear dispersion relations will be given. In the main part of the boundary layer, region 𝑅 1 , we define the coordinate 𝑦 = 𝜖 6 𝑌 , where 𝑌 = 𝑂 ( 1 ) , and introduce a small disturbance of size 𝛿 into the basic flow. The expand disturbance quantities are then expanded as 𝑢 = 𝑈 𝐵 + 𝛿 𝑢 0 + 𝜖 𝛿 𝑢 1 + , ( 3 . 1 ) 𝑣 = 𝜖 𝛿 𝑣 0 + 𝜖 2 𝛿 𝑣 1 + , ( 3 . 2 ) 𝑝 = 𝑃 𝐵 + 𝛿 𝜖 𝑝 0 + 𝜖 2 𝛿 𝑝 1 + , ( 3 . 3 ) 𝑇 = 𝑇 𝐵 + 𝛿 𝑇 0 + 𝜖 𝛿 𝑇 1 + , ( 3 . 4 ) 𝐶 = 𝐶 𝐵 + 𝛿 𝐶 0 + 𝜖 𝛿 𝐶 1 + , ( 3 . 5 ) where 𝑢 𝑖 , 𝑣 𝑖 , and so forth are functions of the boundary layer variable 𝑌 and of the scaled streamwise variable 𝑋 , and 𝛿 is the amplitude of the disturbance which is very much smaller than unity so that terms quadratic in 𝛿 are ignored thereby restricting the analysis to linear stability theory. Substituting (3.1)–(3.5) into the governing equations and solving the resulting system of leading-order equations yields the following first-order solutions: 𝑢 0 = 𝐴 0 𝑈 𝐵 𝑌 , 𝑣 0 = 𝛼 0 𝐴 0 𝑋 𝑈 𝐵 , 𝑇 0 = 𝐴 0 𝑇 𝐵 𝑌 , 𝐶 0 = 𝐴 0 𝐶 𝐵 𝑌 , 𝑝 0 = 𝑃 0 + 𝐺 0 𝑡 𝐴 0 𝑇 𝐵 𝛾 1 + 𝐺 0 𝑐 𝐴 0 𝐶 𝐵 𝜇 1 . ( 3 . 6 ) At the next order, 𝑂 ( 𝛿 2 ) , we obtain the solutions 𝑣 1 = 𝛼 0 𝑐 0 𝐴 0 𝑋 + 𝛼 0 𝑈 𝐵 2 𝑌 0 𝑃 0 𝑋 𝑈 2 𝐵 𝑑 𝑌 + 𝛼 0 𝑈 𝐵 𝐺 0 𝑡 𝐴 0 𝑋 𝑌 𝑌 0 𝑇 𝐵 𝛾 1 𝑈 2 𝐵 𝑑 𝑌 𝛼 1 𝐴 0 𝑋 + 𝛼 0 𝑈 𝐵 𝐺 0 𝑐 𝐴 0 𝑋 𝑌 𝑌 0 𝐶 𝐵 𝜇 1 𝑈 2 𝐵 𝑑 𝑌 𝛼 0 𝐴 1 𝑋 𝑈 𝐵 , 𝑝 1 = 𝑃 1 𝛼 2 0 𝐴 0 𝑌 𝑌 0 𝑈 2 𝐵 𝐷 𝑑 𝑌 0 𝐴 0 𝐺 0 𝑡 𝛼 0 𝑌 𝑌 0 𝑇 𝐵 𝑌 𝑈 𝐵 𝐺 𝑑 𝑌 0 𝑐 𝐴 0 𝐷 0 𝛼 0 𝑌 𝑌 0 𝐶 𝐵 𝑌 𝑈 𝐵 𝑑 𝑌 + 𝐴 1 𝐺 0 𝑡 𝑇 𝐵 + 𝐺 0 𝑐 𝐶 𝐵 𝐺 0 𝑡 𝑌 𝑌 1 𝑇 𝐵 𝑌 ( 𝑌 1 0 [ 𝑃 0 + 𝐺 0 𝑡 𝐴 0 𝑇 𝐵 𝛾 1 𝑈 2 𝐵 ] 𝑑 𝑌 1 ) 𝑑 𝑌 𝐺 𝑐 𝑌 𝑌 1 𝐶 𝐵 𝑌 ( 𝑌 1 0 [ 𝑃 0 + 𝐺 0 𝑐 𝐴 0 𝐶 𝐵 𝜇 1 𝑈 2 𝐵 ] 𝑑 𝑌 1 ) 𝑑 𝑌 𝐴 0 𝐺 0 𝑡 𝐺 0 𝑐 { 𝑌 0 𝑇 𝐵 𝑌 ( 𝑌 1 𝑌 0 𝑐 𝐵 𝜇 1 𝑈 2 𝐵 𝑑 𝑌 1 ) 𝑑 𝑌 + 𝑌 0 𝐶 𝐵 𝑌 [ 𝑌 𝑌 0 𝑇 𝐵 𝛾 1 𝑈 2 𝐵 𝑑 𝑌 1 ] 𝑑 𝑌 } , ( 3 . 7 ) where 𝐴 𝑖 = 𝐴 𝑖 ( 𝑋 ) and 𝑃 𝑖 = 𝑃 𝑖 ( 𝑋 ) ( w h e r e , 𝑖 = 0 , 1 ) are unknown functions representing the displacement and the pressure amplitudes. In the results above, we set 𝐴 𝑖 = 𝐴 𝑖 𝑒 𝑖 𝑋 + c . c , 𝑃 𝑖 = 𝑃 𝑖 𝑒 𝑖 𝑋 + c . c (where c . c denotes the complex conjugate). The lower limit of the integrals, 𝑌 0 , is a non-zero constant introduced for convenience, whose value does not alter the eventual results for wave numbers and frequencies.

In region 𝑅 2 , we define 𝑦 = 𝜖 7 𝑌 with 𝑌 = 𝑂 ( 1 ) and the expansions follow from 𝑅 1 : 𝑢 = 𝜆 1 𝜖 𝑌 + 𝜖 2 𝜆 2 𝑌 2 𝑢 + 𝛿 ( 0 ) + 𝜖 𝑢 ( 1 ) , + 𝑣 = 𝛿 𝜖 𝜖 𝑣 ( 0 ) + 𝜖 2 𝑣 ( 1 ) , + 𝑝 = 𝑝 𝐵 + 𝛿 𝜖 𝑝 ( 0 ) + 𝜖 2 𝑝 ( 1 ) , + 𝑇 = 𝛾 1 + 𝜖 𝛾 2 𝑌 + 𝜖 2 𝛾 3 𝑌 2 𝑇 + 𝛿 ( 0 ) + 𝜖 𝑇 ( 1 ) , + 𝐶 = 𝜇 1 + 𝜖 𝜇 2 𝑌 + 𝜖 2 𝜇 3 𝑌 2 𝐶 + 𝛿 ( 0 ) + 𝜖 𝐶 ( 1 ) . + ( 3 . 8 ) Substituting these equations into the governing equations and solving the resulting equations yield the following first-order solutions: 𝑢 ( 0 ) = 𝜆 1 𝐴 0 , 𝑣 ( 0 ) 𝛼 = 0 𝑝 𝑋 ( 0 ) 𝜆 1 𝛼 0 𝐴 0 𝑋 𝜆 1 𝑝 𝜉 , ( 0 ) = 𝑃 ( 0 ) = 𝑐 0 𝜆 1 𝐴 0 , 𝐶 ( 0 ) = 𝜇 2 𝜆 1 2 𝐴 0 𝜆 1 2 𝜉 + 𝑃 ( 0 ) 𝜉 𝑖 𝐷 0 𝛼 0 𝜆 1 1 , 𝑇 ( 0 ) = 𝛾 2 ( 𝐴 0 + 𝑝 ( 0 ) 𝜆 1 2 𝜉 ) , ( 3 . 9 ) where 𝜉 = 𝑌 ( 𝑐 0 / 𝜆 1 ) . At the next order, the velocity and pressure terms are 𝑣 ( 1 ) 1 = 𝜆 1 𝛼 0 𝑃 𝑥 ( 1 ) + 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 𝜆 1 𝛼 0 𝐴 0 𝑥 𝜉 l n | 𝜉 | + 𝜙 ± 𝜆 1 𝛼 0 𝜉 𝐴 1 𝑥 𝛼 0 𝜆 2 𝐴 0 𝑥 𝜉 2 + 2 𝑐 0 𝜆 1 𝜉 | | 𝜉 | | l n + 𝜙 ± 𝑐 2 0 𝜆 2 1 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 𝜆 1 2 𝛼 0 𝑐 0 𝐴 0 𝑥 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 2 𝜆 1 𝐴 0 𝑥 𝛼 0 𝑐 0 + 𝑖 𝐷 0 | | | 𝜉 l n 2 + 𝐷 0 2 𝛼 0 2 𝜆 1 2 | | | 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 𝐷 0 𝜆 1 𝛼 0 𝐴 0 𝑥 𝛼 0 𝑐 0 + 𝑖 𝐷 0 𝜉 t a n 1 𝛼 0 𝜆 1 𝜉 𝐷 0 𝑖 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 2 𝜆 1 𝐷 0 𝛼 0 𝐴 0 𝑥 𝛼 0 𝑐 0 + 𝑖 𝐷 0 𝜉 | | | 𝐷 l n 0 2 + 𝛼 0 2 𝜆 1 2 𝜉 2 𝛼 0 2 𝜆 1 2 𝜉 2 | | | + 𝜙 ± 𝛾 + 𝑖 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 𝜆 1 2 𝐴 0 𝑥 𝛼 0 𝑐 0 + 𝑖 𝐷 0 t a n 1 𝐷 0 𝛼 0 𝜆 1 𝜉 + 𝜙 ± , 𝑝 ( 1 ) = 𝑃 ( 1 ) + 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 𝐴 0 𝛾 𝑌 + 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 𝛼 0 𝑐 0 + 𝑖 𝐷 0 𝐴 0 | | | 𝜉 l n 2 + 𝐷 0 2 𝛼 0 2 𝜆 1 2 | | | 𝑖 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 𝑐 𝛼 0 𝜆 1 𝛼 0 𝑐 0 + 𝑖 𝐷 0 𝐴 0 t a n 1 𝐷 0 𝛼 0 𝜆 1 𝜉 + 𝜙 ± . ( 3 . 1 0 ) The solutions in this region possess both logarithmic and algebraic singularities as 𝜉 0 . These singularities are smoothed out by the introduction of the critical layer consisting of a thin viscous region situated in the neighbourhood of the critical level 𝜉 = 0 where 𝜙 ± and 𝜙 ± 𝑝 are the phase-shift terms introduced to connect the solutions in the normal velocity and pressure, respectively, on either side of the critical layer.

Solutions in the other regions (namely, the wall layer 𝑅 4 and the outer potential-flow layer 𝑅 5 ) follow in a straightforward manner, and the important solutions of the wall layer are given by 𝑣 0 = 𝑖 𝛼 2 0 𝑃 0 𝑚 𝜔 0 1 𝑚 𝑍 𝑒 𝑚 𝑍 , 𝑝 0 = 𝑃 0 ( 𝑋 , 𝑍 ) , ( 3 . 1 1 ) where 𝑚 = ( 𝛼 0 𝑐 0 ) 1 / 2 𝑒 𝑖 𝜋 / 4 , 𝑦 = 𝜖 8 𝑍 , and 𝑍 is an 𝑂 ( 1 ) coordinate.

In region 𝑅 5 , we set 𝑦 = 𝜖 5 ̂ 𝑦 , where ̂ 𝑦 𝑂 ( 1 ) . The leading-order solutions are given by ̂ 𝑢 0 𝑃 = 0 𝑒 𝛼 0 , ̂ 𝑣 ̂ 𝑦 0 𝑃 = 𝑖 0 𝑒 𝛾 0 ̂ 𝑦 , ̂ 𝑝 0 = 𝑃 0 𝑒 𝛾 0 ̂ 𝑦 , ( 3 . 1 2 ) where 𝑃 0 is an unknown function which describes the disturbance pressure at the outer extreme of the boundary layer. At the next order, the important solutions are ̂ 𝑣 1 𝑃 = 𝑖 [ 1 𝜔 0 𝛼 0 𝑃 0 ] 𝑒 𝛼 0 ̂ 𝑦 , ̂ 𝑝 1 = 𝑃 1 𝑒 𝛼 0 ̂ 𝑦 , ( 3 . 1 3 ) where 𝑃 1 is an unknown function which describes the disturbance pressure at the outer extreme of the boundary layer.

4. Linear Neutral Results and Eigenrelations

In this section, we asymptotically match the solutions in their respective overlap regimes. We will be matching the normal velocities and pressures of the same orders in these respective overlapping regions. The first eigenrelation is found by matching the first-order solutions across the entire boundary layer flow regime to be 𝑐 0 𝜆 1 = 𝐺 0 𝑡 𝛾 1 + 𝐺 0 𝑐 𝜇 1 + 𝛼 0 . ( 4 . 1 ) The matching of second-order pressure components between 𝑅 1 (as 𝑌 ) and 𝑅 5 (as ̂ 𝑦 0 ) yields 𝑃 1 = 𝑃 1 𝛼 0 2 𝐴 0 𝐽 0 𝐺 0 𝑡 𝐴 0 𝐽 3 + 𝐷 0 𝛼 0 𝐽 1 𝐺 0 𝑐 𝐴 0 𝐽 4 + 𝐷 0 𝛼 0 𝐽 2 + 𝐺 0 𝑡 𝐴 1 𝑇 𝐵 + 𝐺 0 𝑐 𝐴 1 𝐶 𝐵 𝐴 0 𝐺 0 𝑡 𝐺 0 𝑐 𝐼 1 , ( 4 . 2 ) where 𝑇 𝐵 = l i m 𝑌 𝑇 𝐵 and 𝐶 𝐵 = l i m 𝑌 𝐶 𝐵 . The constants 𝐼 𝑖 𝑠 and 𝐽 𝑖 𝑠 for 𝑖 = 1 , , 4 are defined in the appendix.

Matching the pressure terms across 𝑅 2 (as 𝑌 ) and 𝑅 1 (as 𝑌 0 ) gives 𝑃 1 = 𝑃 ( 1 ) 𝛾 𝑖 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 𝛼 0 𝜆 1 𝛼 0 𝑐 0 + 𝑖 𝐷 0 𝜙 + 𝐴 1 𝐺 0 𝑡 𝑇 0 𝐵 𝐴 1 𝐺 0 𝑐 𝐶 0 𝐵 , ( 4 . 3 ) where 𝑇 0 𝐵 = l i m 𝑌 0 𝑇 𝐵 and 𝐶 0 𝐵 = l i m 𝑌 0 𝐶 𝐵 . Matching the pressure terms across 𝑅 4 (as 𝑌 ) and 𝑅 2 (as 𝑌 0 ) gives 𝑝 1 = 𝑃 ( 1 ) + 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 2 𝛼 0 𝜆 1 𝛼 0 𝑐 0 + 𝑖 𝐷 0 | | | 𝛼 l n 0 2 𝑐 0 2 + 𝐷 0 2 𝛼 0 2 𝜆 1 2 | | | 𝛾 𝑖 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 𝛼 0 𝜆 1 𝐴 0 t a n 1 𝐷 0 𝛼 0 𝑐 0 𝛾 𝑖 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 𝛼 0 𝜆 1 𝛼 0 𝑐 0 + 𝑖 𝐷 0 𝜙 . ( 4 . 4 ) Matching the normal velocity components between 𝑅 1 (as 𝑌 ) and 𝑅 5 (as ̂ 𝑦 0 ) at the second order gives 𝑃 1 𝑃 = 0 𝑐 0 + 𝛼 0 𝐴 1 𝑈 𝐵 𝛼 0 𝑐 0 𝐴 0 𝛼 0 𝑈 𝐵 𝑃 0 𝐼 2 𝛼 0 𝑈 𝐵 𝐺 0 𝑡 𝐴 0 𝐻 1 𝛼 0 𝑈 𝐵 𝐺 0 𝑐 𝐴 0 𝐻 2 + 𝛼 1 𝐴 0 . ( 4 . 5 ) Matching the normal velocity across regions 𝑅 2 and 𝑅 1 yields 𝐵 0 𝐴 0 𝑋 𝛼 0 𝜆 1 𝐴 1 𝑋 + 𝐿 0 𝐺 0 𝑡 𝐴 0 𝑋 + 𝐿 1 𝐺 𝑐 𝐴 0 𝑋 = 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 𝜆 1 𝛼 0 𝐴 0 𝑋 𝜙 + 𝜆 1 𝛼 0 𝐴 1 𝑋 2 𝛼 0 𝜆 2 𝑐 0 𝐴 0 𝑋 𝜆 1 𝜙 + 𝑖 𝛼 0 𝛾 0 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 2 𝜆 1 𝐷 0 𝛼 0 𝑐 0 + 𝑖 𝐷 0 𝐴 0 𝑋 𝜙 + , ( 4 . 6 ) where the constants 𝐵 0 and 𝐿 𝑖 𝑠 for 𝑖 = 0 , 1 are defined in the appendix. Lastly, a matching of the normal velocity between 𝑅 2 and 𝑅 4 gives 𝑖 𝛼 0 𝑃 0 𝑚 𝑐 0 𝛼 = 0 𝑃 𝑋 ( 1 ) 𝜆 1 + 𝐵 1 + 𝐸 0 𝜆 1 + 𝜔 1 𝛼 1 𝑐 0 𝐴 0 𝑋 + 𝛼 0 𝑐 0 𝐴 1 𝑋 2 𝛼 0 𝜆 2 𝑐 0 2 𝜆 1 𝐴 0 𝑋 𝜙 𝑐 0 𝛼 0 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 𝜆 1 2 𝐴 0 𝑋 𝜙 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 2 𝜆 1 2 𝐴 0 𝑋 𝛼 0 𝑐 0 + 𝑖 𝐷 0 | | | 𝑐 l n 0 2 𝜆 1 2 + 𝐷 0 2 𝛼 0 2 𝜆 1 2 | | | 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 𝐷 0 𝜆 1 2 𝑐 0 𝛼 0 𝐴 0 𝑋 𝛼 0 𝑐 0 + 𝑖 𝐷 0 t a n 1 𝛼 0 𝑐 0 𝐷 0 × 𝑖 𝛼 0 𝑐 0 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 2 𝜆 1 2 𝐷 0 𝛼 0 𝑐 0 + 𝑖 𝐷 0 | | | 𝐷 l n 0 2 + 𝛼 0 2 𝑐 0 2 𝛼 2 𝑐 0 2 | | | + 𝑖 𝛼 0 𝑐 0 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 2 𝜆 1 2 𝐷 0 𝐴 0 𝑋 𝛼 0 𝑐 0 + 𝑖 𝐷 0 𝜙 𝑖 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 𝜆 1 2 𝐴 0 𝑋 𝛼 0 𝑐 0 + 𝑖 𝐷 0 t a n 1 𝐷 0 𝛼 0 𝑐 0 + 𝑖 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 𝑐 𝜆 1 2 𝐴 0 𝑋 𝛼 0 𝑐 0 + 𝑖 𝐷 0 𝜙 ( 4 . 7 ) The relations (4.1)–(4.7) above may be used to eliminate 𝐴 1 , 𝑃 ( 1 ) , 𝑃 0 , 𝑃 0 to obtain a relation which determines the higher-harmonic components of 𝐴 1 . If we restrict our attention to the 𝑒 𝑖 𝑋 components, then, after some algebra, (4.1)–(4.7) lead to 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 𝜆 1 𝛼 0 𝑐 0 𝐴 0 𝑋 𝜙 + 𝜙 2 𝑐 0 2 𝛼 0 𝜆 2 𝜆 1 𝐴 0 𝑋 ( 𝜙 + 𝜙 𝑖 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 ) 2 𝜆 1 𝐷 0 𝐴 0 𝑥 𝛼 0 𝑐 0 𝛼 0 𝑐 0 + 𝑖 𝐷 0 𝜙 + 𝜙 𝑖 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 ) 𝜆 1 𝐴 0 𝑋 𝛼 0 𝑐 0 + 𝑖 𝐷 0 𝜙 + 𝜙 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 ) 2 𝜆 1 𝐴 0 𝑋 𝛼 0 𝑐 0 + 𝑖 𝐷 0 | | | 𝑐 l n 0 2 𝛼 0 2 + 𝐷 0 2 𝛼 0 2 𝜆 1 2 | | | 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 ) 𝐷 0 𝜆 1 𝐴 0 𝑋 𝛼 0 𝑐 0 𝛼 0 𝑐 0 + 𝑖 𝐷 0 t a n 1 𝛼 0 𝑐 0 𝐷 0 + 𝑖 𝛼 0 𝑐 0 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 ) 2 𝜆 1 𝐷 0 𝛼 0 𝑐 0 + 𝑖 𝐷 0 | | | 𝐷 l n 0 2 + 𝛼 0 2 𝑐 2 0 𝛼 0 2 𝑐 2 0 | | | 𝐴 0 𝑥 𝑖 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 ) 𝜆 1 𝐴 0 𝑋 𝛼 0 𝑐 0 + 𝑖 𝐷 0 t a n 1 𝛼 0 𝑐 0 𝐷 0 𝐵 3 𝐴 0 𝑋 + 𝐵 4 𝐴 1 𝑋 = 0 , ( 4 . 8 ) where 𝐵 3 and 𝐵 4 are defined in the appendix. The results for linear theory are now well known and are derived by taking the jump across the critical layer, 𝜙 to be equal to 𝑖 𝜋 . Taking the real parts of equation (4.8) then gives 𝛼 0 𝜆 2 1 2 𝑚 2 𝛼 0 𝜆 2 𝑐 0 2 𝜆 1 𝜋 + 2 𝛼 0 𝑐 0 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 𝜋 + 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 𝜆 1 𝐷 0 𝛾 𝜋 + 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 2 𝜆 1 | | | 𝑐 l n 0 2 𝛼 0 2 + 𝐷 0 2 𝛼 0 2 𝜆 1 2 | | | 𝛼 2 0 𝑐 0 2 𝛾 2 𝐺 0 𝑡 + 𝜇 2 𝐺 0 𝑐 2 𝜆 1 𝐷 0 | | | 𝑐 l n 0 2 𝛼 0 2 + 𝐷 0 2 𝛼 0 2 𝑐 0 2 | | | = 0 , ( 4 . 9 ) where 𝑚 = 𝛼 0 𝑐 0 . However, (4.1) and (4.9) are the crucial eigenvalue relations which fix the neutral wave number to the neutral wave speed.

5. Results and Discussion

To obtain a clearer understanding of the effects of fluid buoyancy and the chemical reaction on the linear stability of the two-dimensional disturbance wave modes, we consider a number of limiting cases when the fluid buoyancy parameters are either large and small and the Damkohler numbers are small. This allows for a detailed examination of the linear eigenrelations (4.1) and (4.9).

We first investigate the limiting behaviour of the neutral eigenrelations as the buoyancy parameters ( 𝐺 0 𝑐 , 𝐺 0 𝑡 ) + . This limiting case corresponds to the increase in the buoyancy force due to an increase in the density difference caused by temperature and chemical concentration differences. Solving the eigenrelations (4.1) and (4.9), we obtain, in the limit ( 𝐺 0 𝑐 , 𝐺 0 𝑡 ) + with 𝐷 0 𝑂 ( 1 ) , 𝛼 0 = 𝛾 2 𝜆 1 𝜆 2 𝛾 1 𝐺 0 𝑡 + 𝜇 2 𝜆 2 𝜆 1 𝜇 1 𝐺 0 𝑐 𝑐 + , 0 = 𝛾 2 𝜆 2 𝐺 0 𝑡 + 𝜇 2 𝜆 2 𝐺 0 𝑐 + . ( 5 . 1 ) These results agree with those obtained by Motsa et al. [7] for the case 𝐺 0 𝑐 = 0 and with those of Shateyi et al. [6] when 𝐺 0 𝑡 = 0 and help to quantify the effects of the combined buoyancy on the normal modes. The results show that, depending on whether the buoyancy terms reinforce or cancel out (e.g., in the case of an endothermic reaction), the normal modes may grow without limit thus hastening the transition to turbulence. A viable transition delay mechanism would be to ensure that the buoyancy terms act contrary to each other so as to reduce the growth of the disturbances.

Solving the eigenrelations (4.1) and (4.9) in the limit 𝐷 0 + with