On Mixed Convection in a Horizontal Channel, Viscous Dissipation and Flow Duality

: The conditions deﬁning a stationary ﬂuid ﬂow may lead to a multiplicity of solutions. This circumstance is widely documented in the literature when mixed convection in a vertical channel or duct is accompanied by an important effect of viscous dissipation. Usually, there are double stationary solutions with a parallel velocity ﬁeld which satisfy given temperature boundary conditions and with a prescribed mass ﬂow rate. However, in a vertical internal ﬂow, the dual solutions can be determined only numerically as they do not have a closed analytical form. This study shows that, in a horizontal channel, stationary mixed convection with viscous dissipation shows up dual ﬂow branches whose mathematical expressions can be determined analytically. The features of these dual ﬂows are discussed.


Introduction
There exists a wide literature pointing out that the governing equations of fluid flow may allow multiple solutions for a given properly conditioned problem. In fact, it is widely assumed that the stationary and fully developed flow of a fluid in a duct or channel with assigned wall conditions for the temperature and for the velocity (typically no slip) is unique provided that a given mass flow rate or pressure gradient is prescribed. This is the case of some elementary solutions for isothermal flow available in all textbooks of fluid dynamics, such as the Poiseuille flow or the Poiseuille-Couette flow, as well as for less known solutions relative to natural or mixed convection heat transfer, such as the Batchelor's profile [1]. However, this widely assumed behaviour is not a general feature of stationary fully-developed internal flows. Significant counterexamples exist when the nonlinear nature of the governing equations for fluid flow emerges in the stationary fullydeveloped regime. There are typical conditions leading to such a situation; one is the non-negligible viscous dissipation effect in a buoyant flow. There are well-documented case studies of this type addressing situations where steady natural or mixed convection dual flows arise in vertical channels or ducts [2][3][4][5][6][7]. A general feature of such solutions is that they are not expressed in a closed analytical form, but they are obtained by numerical methods, usually the shooting method implemented, for example via a Runge-Kutta solver or a power-series technique yielding a recursive determination of the series coefficients.
The analysis of viscous dissipation effects in forced convection and mixed convection flows has a significant impact for both engineering and geophysical applications. Engineering applications are oriented to the analysis and the design of the thermal treatment of fluids, important to the food processing industry and to the chemical engineering analysis of organic and/or polymeric fluid flows. Geophysical applications are devised as typical of the interplay between buoyancy and viscous dissipation effects, as pointed out in the classical papers by Gebhart [8] and by Gebhart and Mollendorf [9]. These authors introduced the dissipation parameter (now widely known as the Gebhart number) as the key parameter governing the phenomenon.
The aim of this paper is to examine the stationary mixed convection flow with viscous dissipation in a horizontal channel, where the channel boundary walls are adiabatic and the velocity field is parallel. It will be shown that there exist dual solutions for every possible assignment of the mass flow rate, expressed in non-dimensional terms through the Péclet number. The striking feature of this flow problem is that the dual solutions are expressed in a closed form as polynomials in the transverse coordinate. This paper mirrors an analogous study carried out for a fluid saturated porous medium governed by Darcy's law [10]. The novel feature of the investigation presented here is that a Newtonian fluid flow is considered instead of the seepage flow in a porous medium.

Governing Equations
Let us consider a horizontal channel bounded by rigid walls at y = 0, H. The y axis is vertical, while the (x, z) axes are horizontal. As shown in Figure 1, the gravitational acceleration is g = −gê y , whereê y is the unit vector along the y axis. The rigid boundary walls are modelled as perfectly adiabatic. We assume the validity of the Boussinesq approximation and we assume the effect of viscous dissipation as non-negligible. Thus, a dimensionless formulation of the local mass, momentum and energy balance equations for a Newtonian fluid flowing in the channel can be expressed as [8,9] ∇ · u = 0, 1 Pr In Equation (3), Φ denotes the dissipation function, where u i and x i are the ith Cartesian components of the velocity vector and of the position vector, respectively, while repeated indices i and j are implicitly summed over.
The dimensionless boundary conditions are formulated as The scaling adopted for the dimensional quantities in order to obtain Equations (1)-(3) and (5) is given by where the subscript "nd" denotes the dimensionless quantities. In Equations (1)-(3) and in the forthcoming equations, such a subscript is omitted for the sake of simplicity. Most of our study, starting from Equation (8), involves dimensionless quantities and, when some dimensional field or quantity is used, this will be explicitly stated to avoid any confusion. In Equation (6), ∆T = αν/(gβH 3 ) and T 0 is the reference temperature employed to evaluate the fluid properties within the Boussinesq approximation. The Prandtl and Gebhart numbers are defined as In Equations (1)- (7), u is the velocity, p is the dynamic pressure, i.e., the local difference between the pressure and the hydrostatic pressure, T is the temperature, and t is time. Furthermore, α is the thermal diffusivity, ν is the kinematic viscosity, β is the thermal expansion coefficient, ρ is the constant reference density and c is the specific heat.
We said that the scaling defined by Equation (6) involves the reference temperature T 0 . The reliability of the model based on the Boussinesq approximation greatly depends on the choice of T 0 . The optimal choice is one that minimises the error in the substitution of the variable density with its reference value. Such a choice of T 0 is the average value over the spatial domain where the flow is studied [11]. The dimensionless spatial region where the flow solution is to be determined can be defined, in terms of (x, y, z), as the to be satisfied by the dimensionless temperature T.

Flow Duality
Equations (1)-(3) and (5) admit stationary two-dimensional solutions, where the velocity field u = (u, v, 0), the dynamic pressure field p and the temperature field T depend only on the coordinates (x, y). Among such solutions, we seek those satisfying the requirement where F(y) is a function satisfying an integral normalisation condition, so that the Péclet number, Pe, expresses the dimensionless average velocity along the x direction. In other words, if we denote by U 0 the dimensional average velocity in a x = constant cross-section, we can write since α/H is the velocity scale employed in Equation (6).

The Velocity Field and the Temperature Field
Equation (9) identically satisfies Equation (1), while Equations (2) and (3) are drastically simplified to where, on account of Equations (4) and (9), the dissipation function Φ is given by We now derive Equation (12) with respect to y, we derive Equation (13) with respect to x and subtract the two resulting equations so that the mixed second-order partial derivatives of p cancel out. Thus, we obtain On account of Equation (9), the boundary conditions (5) yield Let us derive Equation (14) with respect to x by taking into account that, from Equation (15), Φ depends only on y. Thus, we have Equation (16) implies that ∂T/∂x depends only on y, so that both ∂ 2 T/∂x 2 and ∂ 3 T/∂x 3 are identically zero. Therefore, the substitution of Equation (16) into Equation (18) yields Thus, we can conclude that F(y) is a fourth degree polynomial in y.
where c 0 , c 1 , c 2 , c 3 and c 4 are constants. By using Equations (10) and (17), we obtain The general solution of Equation (16) can be written as where G(y) is a function to be determined, and Equation (20) has been employed. On account of Equation (22), the boundary conditions (5) yield which can be satisfied for every x if and only if Then, Equation (21) yields We can rewrite Equations (20) and (22) as The integral constraint (8) is satisfied, when T is expressed by Equation (26), if G(y) has a zero average value, By substituting Equations (15), (16) and (26) into Equation (14), we obtain By employing the expression of F(y) given by Equation (26), we can reformulate Equation (28) as An integration of Equation (29) yields where the condition dG(y)/dy y=0 = 0, prescribed by Equation (24), is identically satisfied. The second boundary condition (24), dG(y)/dy y=1 = 0, is satisfied if the constant c 1 is a solution of There exist two real solutions of Equation (31), if and only if Ge ≤ √ 15, namely When Ge = √ 15, the two solutions coincide as A − = A + = 2 √ 15. By further integrating Equation (30), one can determine G(y) as where the arbitrary integration constant has been chosen so that the integral constraint (27) is satisfied.

The Dynamic Pressure Field
After the determination of the velocity and temperature fields, we can determine the dynamic pressure gradient by employing Equations (12) and (13). In fact, one obtains where Equation (26) has also been employed.

Some Elementary Physics Behind the Dual Flows
One may wonder how the presence of the buoyancy force, together with a nonuniform temperature profile, can be compatible with the absence of a vertical component of the velocity vector, i.e., the absence of an upward motion of the fluid, as declared in Equation (9). The simplest example of a non-uniform temperature profile where the buoyancy force is unable to generate a stationary motion with a nonzero vertical component of the velocity is the basic state of a Rayleigh-Bénard system. We consider a horizontal fluid layer where the plane boundaries are kept at different temperatures with heating from below. This arrangement yields a basic stationary state where the fluid is at rest (no motion either vertical or horizontal) and the temperature field is non-uniform in the vertical direction, with a uniform temperature gradient. Why is the buoyancy force unable to produce an upward motion? The reason is that the buoyancy force is not sufficiently intense to counterbalance the friction force caused by the fluid viscosity if not for small perturbations of the rest state which eventually are damped out in time and do not produce any permanent (stationary) upward flow. Things become different in the Rayleigh-Bénard system only when the temperature difference between the plane boundaries of the fluid layer becomes sufficiently large for the Rayleigh number to exceed its critical value. Then, the rest state may turn into a system of stationary Rayleigh-Bénard cells and a nonzero upward component of the velocity is observed.
The dual flows studied in our paper are stationary and display a zero vertical component of the velocity, despite the nonzero buoyancy force and the non-uniform temperature profile. The physical reason behind this behaviour is exactly the same as that discussed above for the Rayleigh-Bénard system, namely the buoyancy force is not sufficiently intense to counterbalance the friction force caused by the fluid viscosity in a permanent (stationary) manner. As for the Rayleigh-Bénard system, the buoyancy force may eventually overcome the friction force and establish a cellular flow system superposed to the dual parallel flows, but this can happen only when the vertical temperature gradient becomes so intense that the dual flows become thermally unstable. The analysis of the thermal instability is beyond the scope of our study though it will be an interesting development to be pursued in the future. For the physics of the Rayleigh-Bénard system, we mention the review paper by Normand et al. [12]. For the emergence of convection cells in a system similar to that examined in our paper, but involving a fluid saturated porous channel governed by Darcy's law, we refer the reader to Barletta and Rees [10].

Discussion of the Results
An analytical solution of the local balance equations for a Newtonian fluid flowing in a horizontal parallel channel bounded by adiabatic walls has been obtained. The characteristic feature of this stationary and two-dimensional solution is that it predicts a duality of flows for a prescribed assignment of the pair (Ge, Pe). It is important to keep in mind that the Gebhart number, Ge, is the only parameter governing the flow system, as is shown by Equations (1)-(3) and (5), while the Péclet number, Pe, serves to fix the mass flow rate associated with the solution. Strictly speaking, the Péclet number plays the role of a scaling parameter for the dual solutions, so that such dual flows degenerate into a single trivial solution when Pe = 0. Indeed, on account of Equations (9), (26) and (34), we find u = 0, T = 0 and ∇p = 0 when Pe = 0. This is what is to be expected physically as the mass flow rate through the channel causes the frictional heating which, in turn, activates the temperature gradient and the buoyancy force. This chain of physical effects breaks down in the absence of a forced flow rate across the channel.
Beyond the scaling effect of the Péclet number, the phenomenon of flow duality is a consequence of a third parameter, c 1 , which depends only on Ge in a double determination given by Equation (32). In other words, the obtained stationary solution features two branches which, hereafter, will be denoted as A − -branch and A + -branch.

Characteristics of the Velocity Profiles
The deep diversity of features between the A − -branch and A + -branch emerges with the analysis of the velocity profiles. As Pe represents an overall scaling factor for the x-component of velocity, the velocity distribution in a x = constant cross-section can be displayed by the function F(y) and by its parametric dependence on the Gebhart number. This is made evident by Figure 2 where the dual branches A − and A + are illustrated for different values of Ge. A very wide range of Gebhart numbers is accounted for in Figure 2, from the limiting case Ge → 0 to Ge = 2. The latter value of Ge is really huge for most real-world applications except for geophysical applications where the channel width can be sufficiently large to achieve these Gebhart numbers. We bear in mind that Turcotte et al. [13] devised a range of Gebhart numbers up to Ge = 3. Figure 2 features a left-hand frame with the A − -branch velocity profiles and a righthand frame with the A + -branch velocity profiles. The comparison reveals that the two branches yield quite dissimilar profiles for a given Gebhart number. The A − -branch shows slight departures from the Poiseuille velocity profile, which is exactly approached in the limit Ge → 0. This is rigorously proved by considering the lowest orders of a power series expansion of F(y) with respect to Ge. In fact, from Equations (26) and (32) with c 1 = 6 + A − , one has where the zero-order term, and hence the limit of F(y) for Ge → 0, is the Poiseuille profile. On account of Figure 2, the velocity profiles for the A − -branch display a gradual, yet slight, loss of symmetry as Ge increases above zero. The cause of the asymmetry is the buoyancy effect produced by the viscous dissipation. Even when Ge moves close to its largest value compatible with the existence of the dual solutions, Ge = √ 15, the departure from the Poiseuille profile becomes marked but not extremely large, as illustrated in Figure 3. The black line in Figure 3 displays the merging A − and A + profiles for Ge = √ 15. We note that the merging velocity profiles display a bidirectional character with a flow reversal region close to y = 1.  Flow reversal is a permanent property of the A + -branch, as shown by both Figures 2 and 3. On the other hand, the A − -branch displays flow reversal close to the upper boundary only in the narrow range 15/4 < Ge ≤ √ 15, as it can be easily inferred from Equations (26) and (32) with c 1 = 6 + A − . Figure 2 shows that the velocity profiles of the A + -branch are bidirectional. The choice of plotting Ge F(y) in this case is a consequence of the singular behaviour of the A + profile F(y) when Ge → 0. This feature is immediately gathered on evaluating, from Equations (26) and (32), the lowest orders of a power series expansion of F(y) with respect to Ge with c 1 = 6 + A + , By employing Equation (36), one obtains lim Ge→0 Ge F(y) = 60 y (1 − y) (1 − 2 y).
The right-hand side of Equation (37) is a cubic polynomial profile describing an odd function of the shifted coordinate y * = y − 1/2.

Characteristics of the Temperature Profiles
If function F(y) represents the velocity distribution in a x = constant cross-section, the temperature profile at x = constant can be displayed by plotting function G(y) as given by Equations (32) and (33). Figure 4 shows that the A − -branch yields an asymptotic temperature profile in the limit Ge → 0. This temperature profile is that corresponding to the forced convection regime where the Poiseuille profile describes velocity. In fact, this conclusion is supported by Equations (32) and (33) when we evaluate the lowest orders of a power series expansion of G(y) with respect to Ge with c 1 = 6 + A − , G(y) ≈ −3 6 y 4 − 12 y 3 + 6 y 2 − 1 5 Ge The black line in the left-hand frame of Figure 4 displays the fourth-order polynomial that appears in Equation (38) multiplied by Ge. Hence, the choice of plotting the scaled function Ge −1 G(y) for the A − -branch profiles in Figure 4 is meant to capture the forced convection asymptotic case, where the buoyancy force is neglected with respect to the dynamic pressure gradient contribution in the local momentum balance. The A − temperature profiles displayed in Figure 4 for nonzero Ge display a departure from the symmetric temperature distribution in the forced convection limit, Ge → 0. As for the velocity profiles, the asymmetry is induced by the effects of the buoyancy force.
The singular behaviour of the A + -branch for Ge → 0 is highlighted by Equation (39). In the limit Ge → 0, the velocity diverges with Ge −1 , as shown by Equation (36), while the temperature diverges with Ge −2 . This limiting behaviour is the reason for the scaling, Ge 2 G(y), used in the right-hand frame of Figure 4.    (26), show that ∂T/∂y is negative in the lower part of the channel for the A − -branch. Physically, this means a possibly unstable thermal stratification of the fluid that may give rise to a thermal instability of the flow. As already pointed out in Section 3.3, a stability analysis is beyond the scope of this paper, but the possible onset of the instability may be an interesting chance for future developments of this research. If the possibility of a thermally unstable behaviour of the A − -branch is realistic, this phenomenon seems much less conceivable for the A + -branch. In fact, for the A + -branch, an unstable thermal stratification in the lower part of the channel is present, but it is counterbalanced, to a large extent, by a markedly stable thermal stratification in the largest part of the channel. If the A + -branch is the physically unexpected feature of the solution obtained in Section 3, such an unexpected result is unlikely to be ruled out by a stability analysis. A similar conclusion is drawn in Barletta and Rees [10] with reference to the fluid flow in a saturated porous channel.

Are the Dual Flows Compatible with the Boussinesq Approximation?
The governing Equation (2) represents the widely employed Boussinesq approximation of the local mass, momentum and energy balance equations. As pointed out by several authors [14][15][16][17], the Boussinesq approximation is a correction to a fully-incompressible, constant properties flow model where the small effect taken into account is the buoyancy force. As such, the Boussinesq approximation corresponds to an asymptotic regime where a dimensionless parameter, tends to zero. Here, ∆T e f f is the actual maximum dimensional temperature difference occurring in the flow domain. We mention that the conceptualisation of buoyancy-induced flows where the approximate scheme is obtained in the asymptotic case → 0 is already quite clearly laid out in the pioneering paper by Boussinesq [18]. The often overlooked situation is when the solution of the Boussinesq approximate equations is incompatible with the smallness of ε. In rare cases, such a situation may arise even for perfectly reasonable values of the dimensionless parameter(s) governing the flow. In our case, the governing parameter is the Gebhart number which is assumed to be finite, though small enough. By considering the small-Ge expansion of the dimensionless temperature for the A + -branch, Equation (39), we can approximate the maximum dimensionless temperature difference on a x = constant cross-section as where Equation (26) has been also employed. Thus, Equations (6), (11), (40) and (41) lead to an expression of for the A + -branch at small Gebhart numbers, where Fr is the Froude number based on the dimensional average velocity across the channel, Equation (42) discloses the reason why the A + -branch is incompatible with the asymptotic regime → 0 characteristic of the Boussinesq approximation. In fact, the Gebhart number is usually small and the smaller is Ge the larger is according to Equation (42). Therefore, the only chance to obtain a compatibility between the A + -branch and the core assumption of the Boussinesq approximation, i.e., the limit → 0, is considering large Gebhart numbers which, as suggested by Figures 4 and 5, yield definitely smaller differences T max − T min . However, such large Gebhart numbers can be observed only with extremely wide flow systems such as those typical of geophysical or astrophysical applications.

Conclusions
The analysis of stationary two-dimensional mixed convection in a horizontal plane channel bounded by adiabatic walls was performed by taking into account the effect of the viscous dissipation caused by the forced flow rate across the channel. The velocity field was assumed to be parallel and directed along the horizontal x-axis. Under such conditions, the buoyancy force yields an effect exploited as a flow duality: two branches of possible flows exist denoted as the A − -branch and A + -branch.
Dual flows, A − and A + , are allowed for every prescribed pair of input parameters (Ge, Pe), where Ge is the Gebhart number and Pe is the Péclet number, provided that Ge ≤ √ 15. The dual flows were obtained as completely analytical solutions of the local balance equations under the Boussinesq approximation. Such solutions have been expressed as polynomials with respect to the vertical y coordinate.
The dual flows exhibit quite diverse features, with the A − -branch showing up small departures from the forced convection limit of Poiseuille flow at small Gebhart numbers, and the A + -branch displaying bidirectional flows with a singular behaviour as Ge tends to zero.
The compatibility between the dual flows and the Boussinesq approximation was discussed. The proposed arguments led us to a scenario where the dual flows yield a conceivable behaviour, for the A − -branch, which is characterised by small asymmetric departures from the Poiseuille velocity profile. The other branch, A + , has characteristics to a large extent incompatible with the asymptotic nature of the Boussinesq approximation intended as the limiting case where → 0, where = β ∆T e f f with β the thermal expansion coefficient and ∆T e f f the maximum dimensional temperature difference within the flow domain.
A challenge for a future development of this analysis is the investigation of a scheme different from the Boussinesq approximation, which may lead to a self-consistent framework of flow duality, where the features of both branches can be determined in a reliable way. Among the many possibilities, one can employ the Gay-Lussac approximation outlined by Mayeli and Sheard [17] or a more complex fully compressible scheme along the lines drawn, for instance in the analysis carried out by Chenoweth and Paolucci [19].

Conflicts of Interest:
The authors declare no conflict of interest.