FLUIDS 3 , 113902 ( 2018 ) Stability of the Blasius boundary layer over a heated plate in a temperature-dependent viscosity flow

We present complementary numerical and asymptotic studies of the flow over a heated, semi-infinite flat plate for a fluid with temperature-dependent viscosity. Liquid-type viscosities are found to entrain both the velocity and temperature profiles closer to the plate with increasing temperature sensitivity; gas-type viscosities are found to exhibit the reverse effect. A linear stability analysis is presented and we find that increasing the temperature dependence of the fluid (from gasto liquid-type behavior) results in an increased critical Reynolds number to a point of maximum stability. Using an energybalance approach, we determine that this behavior is primarily driven by the inviscid instability of the modified steady flow, rather than being a result of modified viscous instability effects. Application and extension of the results are considered in the context of chemical vapor deposition.


I. INTRODUCTION
The Blasius boundary layer has been a staple of fluid mechanics since its inception in the early 20th century [1].The simple geometry of the flat plate and the relative simplicity of the similarity solutions describing flows over it allows for the theory to be applied to a wide range of practices.As a result, there is still a great output of literature concerning specific types of flows over a flat plate and their effect on the Blasius boundary layer.
The effects of small disturbances on the laminar flow over a flat plate were first investigated by Tollmien [2] and Schlichting [3], who considered the amplification of linear, wavelike disturbances as the primary mechanism for flow instability.A particularly useful result from this work was the ability to identify a critical Reynolds number: the specific point at which the beginning of transition from the laminar to turbulent regime occurs.The existence of Tollmien-Schlichting (T-S) waves were verified by Schubauer [4] after observing the waves for the first time in his heavily damped wind tunnel experiments.
Since then, improved computational and experimental methods, as well as different formulations, have led to considerable variation in the critical Reynolds number across the literature.Jordinson [5] derived a strictly parallel-flow model, presenting solutions to the Orr-Sommerfeld equation that resulted in a critical Reynolds number of 520.A significant discrepancy between this result and experimental data was highlighted by Ross et al. [6], whose experimental study produced a critical Reynolds number of approximately 400.
In an effort to improve agreement with experimental data, many studies have considered nonparallel effects due to boundary-layer growth at the leading edge.Smith [7] utilized triple deck theory to study the lower branch structure of the neutral stability curve.Here, the parallel results appear as the leading terms and nonparallel corrections are produced with successive asymptotic expansions.For high Reynolds numbers, the solutions showed improved agreement with the experimental results of Ross et al. [6] over the purely parallel solutions.Gaster [8] used a method of successive approximations to obtain corrections to the leading order parallel-flow solutions, yielding a critical Reynolds number of approximately 480.Fasel and Konzelmann [9] present direct numerical solutions to a fully nonparallel problem, citing the inconsistency across studies in deeming which nonparallel terms are significant and which are not.The neutral curve data produced by Fasel and Konzelmann [9] is consistent with that of Gaster [8], although in this study emphasis is placed on the importance of using solutions to the Orr-Sommerfeld equation as a first approximation to extend to nonparallel theory.
While linear stability analyses typically aim to model the early stages of transitional flows, numerous studies have highlighted the importance of nonlinearity and non-normality in the prediction of transition to turbulence, as well as their role in the formation of primary instabilities.The e N analysis developed by van Ingen [10] is often used to consolidate nonlinearity within the linear framework, utilizing the linear amplification rates to find empirical agreement with transitional Reynolds numbers observed in experiments.Bertolotti et al. [11] utilize parabolic stability equations to examine both nonparallel and nonlinear effects on the Blasius boundary layer stability, concluding that while both are destabilizing with reference to the linear results, neither is significant enough to account for the discrepancies between neutral curve data obtained from theory and experiments.
Non-normal studies allow for the investigation of global instabilities in the Blasius boundary layer and their role in transition.In response to the assertion of Huerre and Monkewitz [12] that a linear analysis is sufficient for investigating global behavior, Chomaz [13] identifies that beyond the threshold for global instability, the instability behavior becomes rapidly nonlinear.Utilizing front propagation theory and direct numerical simulation (DNS), Cossu et al. [14] examine the nonlinear behavior of disturbances beyond the threshold for global instability, finding that for large-amplitude impulses the globally unstable behavior of the Blasius boundary layer is decisively nonlinear, yet the wave front speed is identical to that of its linear counterpart.As a result, conducting local linear analyses is still of importance when extending to a global stability analysis to inform the nonlinear behavior.Åkervik et al. [15] examines the stability of the Blasius boundary layer through a global eigenmode analysis.Here, they note the non-normal Orr mechanism that triggers T-S type instabilities, showing that this occurs through extraction of energy from the mean flow.This particular result holds interesting relevance to results obtained in this study in Sec.IV, which we discuss in Sec.V.
The study we present here can be considered as the first step in analyzing the stability of two-dimensional temperature-dependent viscosity flows related to the process of chemical vapor deposition; the nonparallel, nonlinear, and non-normal approaches utilized in the aforementioned studies are clearly possible extensions of our work.Though beyond the scope of this paper, the interested reader is directed to the review of non-normal methods and nonlinear stability behavior by Chomaz [13].
The significance of viscous effects on flow instability has led to interest in variable viscosity flows.Strong temperature gradients can cause the viscosity of a fluid to vary immensely, and, as a result, there is a great deal of literature studying temperature-dependent viscosity flows.Kafoussias and Williams [16] study the free-forced vertical flow over a flat plate with an inverse-linear variation of viscosity with temperature.Here, it is found that fluids that become more viscous with increasing temperature create both a narrower temperature profile and velocity profile, while the effect is reversed for fluid viscosities that decrease with temperature.Wall and Wilson [17] consider the effects of two different exponential-type viscosity-temperature relationships on Blasius-flow stability.In both models they note that for an isothermally heated plate, an exponentially decaying viscosity destabilizes the flow, while an exponentially-increasing viscosity is stabilizing.The inverse linear temperature dependence of viscosity is also considered by Jasmine and Gajjar [18] in their study of flow stability over a rotating disk.There they observe the same effects as Kafoussias and Williams [16] for the disk boundary layer, as well as finding that a viscosity that decreases more rapidly with increasing temperature destabilizes the flow.Again, this broadly agrees with the findings of Wall and Wilson [17].In this study, we also consider the inverse-linear relationship between viscosity and temperature-utilized by Kafoussias and Williams [16] and Jasmine and Gajjar [18], among others-established by Ling and Dybbs [19] as accurately modeling the viscosity-temperature behavior of a range of real fluids.
A major motivation for this work is its relevance to chemical vapor deposition (CVD): a microfabrication process where a gas mixture is pumped into a reactor and flows over a heated reactant surface to chemically deposit a thin film of material (see Fig. 1).It is important that laminar flow is maintained in a CVD reactor to promote regular and cohesive film growth, where turbulent or unstable flows could disrupt the deposition process.Longitudinal roll instabilities are well established in buoyancy-driven flows with a free convection element [20], and have been observed in CVD reactors at Rayleigh numbers beyond 1708 (where the Rayleigh number is the product of the Grashof and Prandtl numbers, i.e., the ratio of buoyancy and viscous forces multiplied by the ratio of momentum and thermal diffusivities) [21].From a numerical perspective, while there is a great deal of literature regarding solutions to the laminar flow in a CVD reactor, there is, to the authors' knowledge, no such study on the existence of unstable or transitional flows due to forced convection.Convective instabilities arising from forced flow are perhaps not considered due to a dichotomic view of the flow being laminar or turbulent, and typically the operating Reynolds number of a CVD reactor is in the range of 1 to 200 [21,22], which most would consider to be firmly in the laminar regime.However, modifications to the stability solutions of Blasius-type flows have shown that destabilizing effects are possible, as demonstrated by Wall and Wilson [17], where a critical Reynolds number of approximately 220 is achieved at one extreme of their viscosity temperature-dependence parameter.Hence, it is possible that certain destabilizing criteria may be sufficiently influential to produce primary instabilities within the operational parameters of a CVD reactor.It is widely accepted that the steep temperature gradients in a CVD reactor cause variation in all physical properties of the fluid [21][22][23], supporting the need for a temperature-dependent viscosity model.This study serves as the first step towards a larger stability analysis of flows in CVD reactors.
Horizontal CVD reactors often feature a shallow susceptor incline to improve growth uniformity [21].Griffiths [24] studies the stability of shear-thinning flows over a wedge of 45 • through a Falkner-Skan-type formulation, finding them to be significantly more stable than Blasius flow solutions.It is possible that unstable flows that exist within the operating range of a CVD reactor could be stabilized by increasing the angle of inclination.It is well established, however, that free-convective longitudinal rolling instabilities can be induced over heated inclined plates [20,25], so the angle of inclination is kept low in CVD reactors to subdue these instabilities.The scope of this paper does not include angle of incline as a parameter, though future work will examine the effect of a small incline on forced convective instabilities, specifically.
The steep temperature gradient within a CVD reactor would suggest that buoyancy effects are significant.In mixed convective flows, a buoyancy parameter often considered is the ratio between the Grashof number and the square of the Reynolds number.The formulation used in this study is consistent with a flow of low Grashof number and a buoyancy parameter much less than 1, such that free convection and buoyancy are considered negligible.Whilst this can be representative for CVD flows under certain conditions, this approach is deemed necessary in order to focus specifically on the forced convective processes under CVD-like conditions.
In Sec.II of this study we present numerical solutions to the laminar flow over a horizontal, heated flat plate for both liquid-and gas-type temperature-viscosity relationships.We then consider the effect of a small disturbance on the velocity, pressure and temperature profiles in Sec.III, and develop this into comprehensive asymptotic and numerical linear stability analyses.In addition to this an integral energy analysis in presented in Sec.IV with the aim of providing physical insights into our results.In Sec.V we discuss the physical interpretation of these results, as well as briefly considering their impact when applied to CVD and the expansion of this work into a detailed model of flow stability in a horizontal CVD reactor.

A. Formulation
We consider a steady, incompressible, Newtonian fluid flowing over a semi-infinite flat plate with velocity U * = (u * , v * ).Here u * and v * are the velocities in the streamwise and wall-normal directions x * and y * , respectively.The temperature of the fluid is denoted by T * , and the plate is heated to a fixed temperature T * w .Note that throughout our analysis asterisks indicate dimensional quantities.The system is governed by the following mass, momentum, and heat continuity equations: Here, ρ * is the fluid density, t * is time, p * is pressure, C * p is the specific heat capacity of the fluid, and k * is the heat diffusion constant.The viscous stress tensor is given by τ * = μ * γ * , where γ * = ∇ * U * + (∇ * U * ) T is the rate-of-strain tensor.The heat continuity equation (1c) is coupled to (1a) and (1b) via the temperature dependence of the fluid viscosity, , where μ * ∞ and T * ∞ represent the free-stream viscosity and temperature, respectively.The characteristic constant of the fluid ε * represents the sensitivity of viscosity to changes in temperature.For a heated plate, i.e., T * w > T * ∞ , as will be considered here, the case of ε * < 0 represents the viscothermal behavior of a gas, while the case when ε * > 0 represents that of a liquid.Note that setting ε * = 0 yields a constant viscosity and uncouples (1c) from (1a) and (1b).
The system is nondimensionalized by introducing the following velocity, length, time, pressure, and temperature scales: where The non-dimensional system of equations becomes ∂u ∂x + ∂v ∂y = 0, (2a) is the Reynolds number, and Pr = C * p μ * ∞ /k * is the Prandtl number.In the proceeding work, we set Pr = 0.72, which is consistent for a range of gases including air.
We now rescale the problem to the boundary-layer region close to the surface of the plate.It is convenient at this point to introduce the small parameter σ = Re −1/8 ; then we have that The boundary-layer equations, at leading order, are then where we have assumed that the pressure is constant across the boundary layer and the subscript notation denotes a partial derivative.
The system (3) is solved subject to the no-slip condition U (Y = 0) = V (Y = 0) = 0, and the far-field free-stream matching condition U (Y → ∞) → 1.In addition to this, due to our choice of nondimensional variables, we have that (Y = 0) = 1 and (Y → ∞) → 0. We define the boundary layer coordinate η = Y/ √ x, and introduce the following similarity variables: where the primes indicate differentiation with respect to η. Substitution of the above into (3) results in the following set of mean flow equations: The boundary-value problem ( 4) is solved using a fourth-order Runge-Kutta integrator alongside a two-dimensional Newton-Raphson shooting routine, subject to the boundary conditions The solutions reached convergence at η = 20 to a tolerance of 10 −8 , and we regard this as numerically representative of free-stream conditions (η → ∞).
We note the effect of ε on the boundary-layer thickness δ * by measuring the Blasius constant δ, given by This definition is also of importance to the analysis presented in Sec.III.

B. Solutions
From Fig. 2(a) we see that the fluid velocity U increases from stationary at the plate surface to the free-stream velocity at some distance η, which varies with ε.As ε is increased, the mean velocity profile is entrained towards the plate surface.This entrainment results in a steeper velocity gradient, as can be seen in Fig. 2(b).From Table I we see that increasing ε also decreases δ.Due to the dependence of δ on U , it can be inferred that the narrowing of the mean velocity profile is akin to the narrowing of the boundary layer.
We note the different shapes of the velocity profiles for positive and negative values of ε.For ε < 0, the velocity gradient is shallow both close to the plate surface and at the boundary-layer edge; while for ε > 0, this behavior is only exhibited at the boundary-layer edge.Figure 2(b) shows this difference in greater detail.For ε 0, the maximum velocity gradient occurs at the plate surface.Conversely, when ε < 0 we observe an inflectional U profile, where the velocity gradient increases to a maximum at some distance away from the plate surface.
When ε = 0, the velocity profile is seen to increase linearly in the region close to the plate, i.e., U is constant.The inflectional U profiles exhibited for the cases when ε < 0 show that this region is moved from the wall towards the center of the boundary layer, while when ε > 0 the profiles do not exhibit any regions where the velocity gradient is constant.As a result, the temperature-independent case, ε = 0, appears to be a limit for linear behavior (U vs η) in the region close to the boundary wall.
Figure 3(a) shows a similar movement towards the plate surface with increasing ε, which we interpret as a narrowing of the thermal boundary layer.This narrowing shows an increase in the temperature profile gradient similar to that seen in Figs.2(a) and 2(b), although the impact of increasing ε is less significant.Figure 3(b) shows the viscosity profiles for a range of temperature dependencies.Variation when ε < 0 results in a much more drastic change in the viscosity profile, and for these values the viscosity decreases rapidly from the plate surface to the free stream.

III. LINEAR STABILITY ANALYSIS
We now assess the stability of the boundary layer when subject to a small disturbance.We seek to consolidate the results of a lower branch asymptotic stability analysis with a numerical approach producing the full neutral stability curves.In flat plate boundary layers the most amplified T-S disturbances appear near to the lower branch of the neutral curve.One has to consider the lower branch structure of the neutral curve in order to investigate the structure of the T-S waves in the near-wall viscous layer.That is our interest here, and our analysis is presented in the following subsection.4. Schematic diagram showing the lower branch structure of the variable viscosity boundary layer flow considered here.The zones I, II, and III denote the upper, main and lower decks respectively.The shaded area indicates the boundary layer region whereas the unshaded area indicates the inviscid region where the base flow matches with that of the free stream.Adapted from Griffiths et al. [26].

A. Asymptotic analysis
In order to conduct the lower-branch asymptotic stability analysis, we assume that the Reynolds number is large so that our parameter σ is small.As in the constant viscosity case (see Smith [7]) we find that the linear disturbances are governed by a triple deck structure on a streamwise length scale of O(σ 3 ).The upper, main, and lower decks are found to be of thickness O(σ 3 ), O(σ 4 ), and O(σ 5 ), respectively (see Fig. 4).The analysis in the upper and main decks is largely similar to that of Smith [7]; however, we find that the viscous lower deck analysis is considerably modified.This is due entirely to the nonconstant, temperature-dependent viscosity function.
In the first instance we make no assumptions regarding the parallel/nonparallel structure of the flow and perturb the base flow as p(x, y, t ) = P 0 + p(x, y, t ), (6c) where P 0 = const.We substitute the above into (2) and retain only the linear perturbation terms to arrive at the governing disturbance equations that are given in Eq. (A1).Due to the scaling of the lower-branch mode, we consider disturbances proportional to where α is the streamwise wave number and β is the frequency of the disturbance.We restrict our attention to neutral disturbances and expand α and β as In the subsequent analysis we adopt a multiple-scales approach whereby the partial x derivative terms of the perturbation quantities in Eq. (A1) are replaced by ∂/∂x + iα/σ 3 .
We introduce the upper deck wall-normal coordinate ȳ = σ −3 y = O(1), and note that here U − 1 = V = = 0.In this deck the we expand the disturbance quantities as In the main deck Y = σ −4 y = O(1), and we find that Substitution of the above into (A1) produces sets of leading order and next order equations in the upper and main decks; for completeness these can be found in the Appendix.Matching between the decks, we find that the solutions at leading order are given by where F 1 = α 1 A 1 e −α 1 ȳ , and A 1 = A 1 (x) is an unknown amplitude function.At the next order we obtain where ], and A 2 = A 2 (x) is another unknown amplitude function.The definitions of the integral functions J X i can also be found in the Appendix.The main deck solutions are then used to match the solutions in the viscous lower deck.We introduce the lower deck wall-normal coordinate Z = σ −5 y = O(1), and note that here the base flow has the form x, and χ = (0)/ √ x.Here we note the emergence of O(σ 2 ) terms in the expansion for U that do not appear in the Newtonian case when ε = 0.
In the lower deck we expand the disturbance quantities as At leading order the equations to solve are then After some manipulation we find that the solutions that match with the those in the main deck are where Ai is the decaying Airy function and The definitions of integral functions I X i can again be found in the Appendix.In order to simplify the analysis in the lower deck, we choose to fix the Prandtl number such that Pr = μ−1 0 = 1 + ε.Solutions are obtainable when Pr = μ0 ; however, this adds significant unnecessary complexity to the problem at the next order.
At the next order the equations to solve are Here we see the emergence of numerous "non-Newtonian" terms owing to both the added complexity of the linear disturbance equations (A1) and the O(σ 2 ) terms in the expansion of the streamwise base flow.In order to determine the second-order eigenrelation, we need only solve for P 2 and U 2 .Matching with the main deck, it is trivial to show that The solution for U 2 is much more involved due largely to the additional terms appearing in Eq. (9b).However, after some lengthy calculations one finds that , where B 2 = B 2 (x) is an unknown amplitude function.We eliminate A 1 (x) from the first-order solutions and obtain the leading order eigenrelation This relation is largely the same as that of Smith [7], differing only by the scaling of the viscosity function at the wall.In order for α 1 to be real (since we are considering neutrally stable modes) we find that ξ 0 ≈ −2.2970i 1/3 = −c 1 i 1/3 and Ai (ξ 0 )/I ∞ 0 ≈ 1.0003i 1/3 = c 2 i 1/3 .Having determined this expression for α 1 , the result for β 1 follows since β 1 = −qα 1 ξ 0 .
After considerable algebra one finds that the second-order eigenrelation is then where It is interesting to see that terms owing from the temperature dependency of the viscosity function appear in the eigenrelation at this order.In an equivalent study of non-Newtonian shear-thinning fluids, Griffiths et al. [26] note that additional viscous terms do not appear in the eigenrelation until the fifth order, the same order at which nonparallel effects are first encountered.This suggests that, at least for the viscosity-temperature relationship considered here, variational viscosity effects will significantly alter the stability properties of the flow when one considers a parallel-flow numerical investigation.
Equating real and imaginary parts of (10), we are able to determine β 2 such that α 2 is real.The singular integral J is computing using Hadamard regularization, and the interested reader is referred to Griffiths et al. [26] for details regarding this calculation.In order to present our results in a manner consistent with the upcoming numerical analysis, we choose to plot the frequency of the disturbance against the Reynolds number.We introduce the frequency parameter where = ε μ0 (0)/2U (0) (we note that this term is identically zero in the constant viscosity case).The expression above represents two terms of the asymptotic expansion of the neutrally stable lower-branch mode.It is clear to see from this result how the stability of the flow at large Reynolds numbers depends on the initial basic structure of the mean flow.In Fig. 5 we plot F against R for a range of ε values.The flow in unstable in the region above the curves.Hence, as the value of ε decreases, we see that the region of instability is significantly increased and the lower-branch mode is destabilized.These results at large Reynolds number suggest that the flow will be marginally stabilized for value of ε > 0 but markedly destabilized when ε < 0. The results presented within this sub-section provide some insight into the stability characteristics of these flows when the Reynolds number is large.In order to be able to make comments regarding the transition from stable to unstable flow at moderate Reynolds numbers we must solve the governing set of linear disturbance equations numerically.Our methodology for doing so is presented in the following sub-section.

B. Numerical analysis
In order to derive the governing set of disturbance equations used in the forthcoming numerical analysis, we nondimensionalize the quantities in Eq. ( 1) with the following velocity, length, time, pressure, and temperature scales: We recall that the choice of scaling is identical to that used in Sec.II with the exception of the length and time scales, which are now scaled in relation to the boundary layer thickness δ * .The perturbed base flow ( 6) is then substituted into the resulting system of equations.We impose a "parallel-flow approximation" whereby the disturbances are considered to occur sufficiently downstream (x 1) such that streamwise boundary-layer growth can be ignored.As a result, the mean flow variables U and are now functions of y alone and we assume that V 1.After linearizing with respect to the perturbation quantities, we arrive at a system of PDEs that are separable in x and t.It is then natural to express the perturbation quantities in terms of normal modes, ĉ = c(y) exp[i(αx − ωt )].Here, α = α r + iα i is the complex streamwise wave number and ω is the real disturbance frequency.
The resulting system of stability equations is then and we note that in the limit as ε → 0 this system can be reduced to the classical fourth-order Orr-Sommerfeld equation.
The system (12) can be solved as a quadratic eigenvalue problem of the form A 2 α 2 + A 1 α + A 0 = Q, where Q = ( ũ, ṽ, p, T ) T is the vector of eigenfunctions, and the quantities A j are matrices containing the coefficients of the O(α j ) terms.The eigenfunctions are then computed according to the boundary conditions which represent the no-slip condition at the wall [the condition for ṽ arises from (12a)], as well as the assumption that fluctuations in temperature relative to the temperature of the wall are negligible.We also assume that all disturbances decay away from the wall.The neutral temporal and spatial stability solutions are obtained via a Chebyshev polynomial discretization method.An exponential map is used to transform the Gauss-Lobatto collocation points into the physical domain.The stability equations are then solved as primitive variables over 100 collocation points distributed between the upper and lower boundaries, with the exception of the conditions described in Eq. ( 13) which are imposed at y = 0 and y = y max , where suitable mean-flow convergence is again found at y max = 20.
Figure 6 shows neutral stability curves (α i = 0), for a range of temperature dependencies, where the curve itself represents solutions to the eigenvalue problem.The region enclosed by the curves are indicative of linearly unstable Reynolds numbers and frequencies.Note that the frequency parameter, F = ω/R, is consistent with the definition provided in Sec.III A. It can be seen that the range of unstable frequencies and Reynolds numbers is reduced with increasing ε.The critical Reynolds number R c , the lowest value of R at which instability is observed, increases with ε.Over the range of ε values studied here, we see that viscosities that increase with temperature (−1 < ε < 0) are less stable than a constant viscosity fluid, while viscosities that decrease with increasing temperature (ε > 0) are more stable.However, this stabilizing effect appears to diminish as ε is increased.
Figure 7 presents a logarithmic comparison between the neutral modes (R vs F ) predicted by the asymptotic and numerical analyses.We observe strong quantitative agreement between the two approaches for both the temperature-independent and temperature-dependent flow cases, particularly at high Reynolds number.In both cases the asymptotic analysis predicts instability at lower frequencies when compared to the numerical results.
The relationship between ε and R c is examined over an extended range of ε in Fig. 8.It can be seen that a point of maximum stability occurs at ε = 1, beyond which the flow begins to destabilize, and increasing temperature sensitivity beyond ε ≈ 4.2 results in a flow that is less stable than one that is temperature independent.
We now examine the behavior of the disturbance eigenfunctions in response to a change in ε.We take this opportunity to introduce a perturbation stream function,  as ε is increased.This may be in response to the narrowing of the corresponding mean flow profiles observed in Figs. 2 and 3. A similar effect was observed by Wall and Wilson [17], who described it as a "compression" effect and attributed it to the mean flow response of their temperature-dependent viscosity model.
It is observed in Fig. 10 that the profile shape close to the disk surface differs for ε < 0. The | ũ| profile maximum flattened at ε −0.25 and −0.5, and is divided into two distinct maxima for ε = −0.75.For ε > 0, the maximum in | ũ| is significantly increased.
In Figure 11, we observe a nonmonotonic relationship between the maximum value of | ṽ| and increasing ε as that seen in Fig. 8.However, this turning point appears in the range 0.25 < ε < 0.75.We speculate on the physical significance of this in Sec.V.
Beyond the narrowing effect also exhibited by the other eigenfunctions, Fig. 12 does not show significant variation of the | T | profile with increasing ε.The impact (or lack thereof) of the temperature eigenfunction on the flow stability is examined in further detail in Sec.IV.

A. Energy balance equation
To better understand the underlying mechanisms of instability, we now consider the energy balance of the system.We multiply streamwise and wall-normal perturbation PDEs by û and v, respectively, and then sum the result to form an equation representing the energy transfer processes of the system.We define the following variables: where ê is the kinetic energy of the two disturbance quantities, q is the disturbance vorticity, and ŝ is defined purely for convenience.
Following some manipulation, we arrive at the resulting energy equation where μ = −ε μ2 T .
Averaging the perturbations over a single time period and integrating up through the boundary layer results in the following integral energy equation: The terms in I represent energy production due to Reynolds stresses (EPRS), terms in II represent energy dissipation due to viscosity (EDV), and terms in III represent additional terms arising from variable viscosity (AVV).As such, when ε = 0, the terms in III vanish.

B. Results
We solve the integral energy integral numerically and observe that the only non-negligible terms (on the left-hand side) are I and the second term of II.It is of particular significance that all terms in III are consistently negligible for the entire range of ε values considered; this shows that the flow stability is not impacted by disturbances to the temperature and the viscosity; rather the stability affects are determined by inviscid responses to changes in the mean-flow profiles.
Given that the perturbations have the normal mode form, the remaining terms are then normalized with respect to ∞ 0 U ẽ + ũ p dy, giving where ẽ = (1/2)( ũ2 + ṽ2 ), q = iα ṽ − ũ , and x ỹ = x ỹ + x ỹ (here indicate complex conjugate).Any eigenmode of a disturbance is amplified when energy production outweighs energy dissipation.The above would then return α i < 0, which is consistent with our definition for linear instability.
Figure 13 shows the change in energy contributions over the range of ε values for R = R c (ε) + 200.The most unstable frequency is used in each case.Solving the energy balance equation at R c + 200 allows the growth rate relative to increasing R to be assessed consistently for each value of ε.Note that the AVV terms have been included to demonstrate that they are consistently negligible.It is observed that when ε < 0 the disturbances are significantly more amplified.The energy produced by the disturbance then decreases as the temperature dependence is increased.Interestingly, the energy profiles exhibit little to no change in behavior beyond the destabilizing turning point outlined in Fig. 8 (i.e., ε > 1).This suggests that although the critical Reynolds number decreases for ε > 1, the growth rates of these disturbances remains near constant.
We observe that the energy loss from EDV decreases for increasing ε.Having established that the influence of ε on the flow stability is via changes in the mean flow, we now manipulate EDV into a constant viscosity component (μ ≡ 1) and a temperature-dependent component: Figure 14 shows how these two components interact over the range of ε values.Although they are notionally viscous dissipation terms, the TDD terms do in fact act as energy production terms when ε > 0, reducing the overall viscous dissipation.We also observe that the ND profile decreases proportionally, leading to the asymptotic behavior of the EDV profile.

V. DISCUSSION AND CONCLUSIONS
We have assessed the onset of instability in the flow of a temperature-dependent fluid over a horizontal flat plate.From solutions to the mean flow in Sec.II, we observe that the boundary layer formed over a horizontal flat plate is broader for gas-type viscosities that increase with temperature, i.e., when ε < 0, than that of a fluid with uniform viscosity, ε = 0.In this case, the increased wall viscosity induces a greater wall shear stress, and the initial deceleration experienced by the fluid near to the plate surface is also greater.The increased viscous forces reduce the flow velocity throughout the boundary layer, causing it to accelerate to that of the free stream more slowly, resulting in a broader boundary layer.
Conversely, setting ε > 0 models a liquid-type viscosity that decreases with increasing temperature.Here there is a reduced wall shear stress, which results in less of a deceleration of the fluid near the plate surface.The viscosity throughout the boundary layer is less than the cases when ε 0, which allows the fluid to accelerate to the free-stream velocity more rapidly, resulting in a narrower boundary layer.
The linear stability analyses in Sec.III show that when ε < 0 the flow is strongly destabilized, while viscosities for which ε > 0 are shown to be stabilizing.There is good agreement between the asymptotic and numerical studies at large Reynolds numbers.
We showed in Sec.IV that the stability of this system is influenced by inviscid responses to the modified mean flow.An interesting parallel is drawn here to the global stability analysis conducted by Åkervik et al. [15], where it is found that the so-called "Orr structures" produced in the Blasius boundary layer gain their energy via the mean shear, which is returned once they decay, triggering T-S type instabilities.This prompts interest in the extension of this work to a global analysis, where the variation in mean flow profiles with temperature dependence seen here could impact the formation and decay of these Orr structures.A comparative study is then also possible, as Åkervik et al. [15] show that local convective stability data is recoverable from specific global eigenmodes.
It is appropriate then to interpret the variation in neutral stability curves seen in Fig. 6 as a result of the changes seen in Figs.2-3.In the case of the broader boundary layer produced by ε < 0, the inflectional U profiles shown in Fig. 2(b) result in greater amplification of the disturbances via the energy production term in Eq. ( 14), as seen in Fig. 13.
The reverse should then be true of a narrower boundary layer, where the velocity profiles increase more uniformly and hence stability is increased.However, a point of particular interest is the appearance of a point of optimum stability for a given temperature dependence, as established in Fig. 8.For mean flow solutions beyond those presented, the Blasius constant continues to decrease, suggesting that the boundary layer continues to narrow and-by the physical interpretation established so far-that this would continue to produce more stable flows, where perhaps R c would behave asymptotically as in other studies [24].Instead, we observe that further narrowing of the boundary layer results in a reduced critical Reynolds number.
The energy production term in Eq. ( 14) is dependent on U .We can see from Fig. 2(b) that increasing ε narrows this profile while increasing the peak at the wall, inducing a steeper shear rate gradient.This change in the U profile must be favorable up to the point of maximum stability, reducing the energy production term and thereby creating a more stable flow.Perhaps beyond this, the differences in velocity through the boundary layer are so great that it renders the flow more susceptible to instability.
We recall from Fig. 11 that the wall-normal eigenfunction also exhibits nonmonotonic behavior with respect to increasing ε, and note that ṽ is also present in the energy production term.It is therefore possible that the increase in the maximum of the | ṽ| profile for ε > 0.5 is a significant mechanism for destabilization of the flow with further increase of ε.Further investigation into the specifics of this point of maximum stability is required.
We briefly consider our results in the context of CVD reactor flows.Recall that the temperature dependence of a gas viscosity is represented by ε < 0, and that the upper limit of operation of a CVD reactor returns flows with R ∼ 200.We see from Fig. 8 that, in the lower limit of ε = −0.75,R c ≈ 70.For context, based on its viscosities at 300 and 1300 K under atmospheric pressure [27], dry air has a value of ε ∼ −0.62 according to our model.It is therefore feasible that a gas mixture used in CVD could have a viscosity-temperature relationship within the parameter range where linear instability appears for R < 200.
There is significant room for development of the current formulation to model flows within a CVD reactor.Jensen et al. [22] identifies the need for temperature-dependent density to reflect gas expansion in the reactor.However, Kays [28] notes that by approximating constant Pr and C p (both being parameters that show little variation with temperature), other temperature-dependent gas properties, e.g., ρ, can be explained as "compensating" for each other.In such a case, the same mean flow equations derived in this study will apply, and the incompressible approximation is sufficient.The extension of this work to form a comparative study that fully incorporates the temperature dependence of each physical property of the fluid is necessary to confirm which of these approximations are justified.
The reaction mechanisms in CVD processes are multistaged and complex, and the most rigorous analysis would require multiple-interacting-species continuity equations for each stage.However, reactants in a CVD reactor are typically very dilute compared to carrier gases, so a "no gas-phase reactions" approximation can be made to allow single diffusion equation to be used to track all chemical species [21].Consistent with this approximation, the Soret effect should also be utilized to highlight the different diffusive responses of the particle mixture to the temperature gradient.ACKNOWLEDGMENTS R.M. would like to thank the EPSRC for funding this research (award reference 1658206).The constructive comments of two anonymous referees are greatly appreciated.

APPENDIX 1. Linear disturbance equations
The governing set of linear disturbance equations are as follows: where V = σ 4 V .Here we note that the fist-order Taylor expansion of the function μ has been used to retain the mean flow form:

10 FIG. 2 .
FIG. 2. (a) Mean flow velocity profiles and (b) velocity gradient profiles for ε = −0.75 to 0.75, with each plot representing an increment of 0.25.The temperature-independent (ε = 0) solution is represented by the black curve.The solutions have been truncated to η = 10.

10 FIG. 3 .
FIG. 3. (a) Mean flow temperature profiles and (b) mean flow viscosity profiles for ε = −0.75 to 0.75, with each plot representing an increment of 0.25.The temperature-independent (ε = 0) solution is represented by the black curve.The solutions have been truncated to η = 10.

FIG. 5 .
FIG.5.Asymptotic results for the neutrally stable lower-branch mode of ε = −0.75 to 0.75, with each line separated by an increment of 0.25.The frequency of the disturbance F is plotted against the Reynolds number R on a log-log scale.

4 FIG. 6 .
FIG.6.Neutral stability curves for ε = −0.75 to 0.75, with each curve separated by an increment of 0.25.Here, F = ω/R is the scaled disturbance frequency, and the axis is scaled logarithmically.The regions enclosed by the curves indicate Reynolds numbers and frequencies for which disturbances create instability in the flow.
as done by Wall and Wilson[17].

3 FIG. 7 .FIG. 8 . 10 FIG. 9 .
FIG. 7. Agreement between the asymptotic theory and the numerical results.The Newtonian solutions are presented in black and the case when ε = −0.28(and hence Pr = 0.72 in the asymptotic framework) is presented in blue.The large Reynolds number asymptotic results are given by the dashed curves.Amplitudes of the perturbation stream function and streamwise, wall-normal, and temperature eigenfunctions for a range of ε values are plotted in Figs.9-12.The eigenfunctions are assessed at R c (ε) + 200.In each case, it can be seen that the profile maximum moves closer to the plate surface

10 FIG. 10 .
FIG.10.Streamwise eigenfunction profiles for ε = −0.75 to 0.75, with each plot representing an increment of 0.25.The temperature-independent (ε = 0) solution is represented by the black curve and all profiles are normalized on the maximum of this profile.The solutions have been truncated to y = 10.

10 FIG. 11 .
FIG.11.Wall-normal eigenfunction profiles for ε = −0.75 to 0.75, with each plot representing an increment of 0.25.The temperature-independent (ε = 0) solution is represented by the black curve and all profiles are normalized on the maximum of this profile.The solutions have been truncated to y = 10.

10 FIG. 12 . 15 FIG. 13 .FIG. 14 .
FIG.12.Temperature eigenfunction profiles for ε = −0.75 to 0.75, with each plot representing an increment of 0.25.The temperature-independent (ε = 0) solution is represented by the black curve and all profiles are normalized on the maximum of this profile.The solutions have been truncated to y = 10.
1) q2 dy TDD , where ND indicates Newtonian dissipation, while TDD indicates temperature-dependent dissipation.The separation of EDV into these separate terms allows for further insight into how the variable viscosity influences flow stability.

TABLE I .
Computed mean flow data for a range of temperature dependencies.