Hexagonal Cell Formation in Darcy–Bénard Convection with Viscous Dissipation and Form Drag

: Recent interest in the effects of viscous dissipation on convective ﬂows in porous media has centred almost exclusively on forced convection ﬂows. In this paper, we investigate the manner in which it affects the onset and early stages of convection in Darcy–Bénard convection. A weakly nonlinear theory is described brieﬂy, and it is shown that hexagonal cells are preferred over rolls when the Rayleigh number is sufﬁciently close to 4 π 2 . At higher Rayleigh numbers, two-dimensional rolls are preferred. When weak form drag is included, then subcritical convection eventually disappears as the Forchheimer parameter increases, yielding a highly novel situation wherein hexagonal convection arises supercritically. The range of stability of hexagons is found to increase.


Introduction
The problem of the onset and development of convection in a uniform horizontal layer heated from below is known as the Darcy-Bénard problem or the Horton-Rogers-Lapwood problem, and it has been studied widely since the pioneering works of Horton and Rogers [1] and Lapwood [2]. Detailed accounts of the development of the subject may be found in the reviews by Nield and Bejan [3], Rees [4], Tyvand [5] and Barletta [6].
In its classical form Darcy's law applies, which means that there are no extensions such as those represented by the boundary and form-drag terms, an anisotropy permeability or temperature-dependent properties. Likewise, the heat transport equation does not include thermal dispersion terms, an anisotropic diffusivity, a temperature-dependent viscosity or local thermal nonequilibrium effects. In this original form, and in an unbounded layer, the critical Rayleigh number is very well-known to be 4π 2 and the associated wavenumber is π [1,2]. Weakly nonlinear theory then shows that two-dimensional rolls are favoured; three-dimensional planforms such as rectangular or hexagonal cells are unstable and therefore unachievable in practice. There are contexts when this straightforward scenario is altered -these include layers with discrete sublayers (square cells; see Rees and Riley [7]), layers with imperfectly conducting boundaries (also square cells; Riahi [8] and Rees and Mojtabi [9]), layers with spatially varying boundary temperatures (rectangles, wavy rolls; see Rees and Riley [10,11]), and layers with internal heating (hexagonal cells; see Tveitereid [12]).
In this paper, we consider the combined effects of viscous dissipation and form-drag (also called the Forchheimer effect) on weakly nonlinear convection in an otherwise classical Darcy-Bénard layer. While one could have chosen a different pair of "effects", these two are found to interact in a somewhat unusual manner, and therefore it is these which will be considered. At present, there is only one paper to our knowledge that considers the effect of viscous dissipation in the nonlinear regime in the context of thermoconvective instability, and that is the analysis by Celli et al. [13] who consider mixed convection in a Darcy-Bénard layer. Their externally-induced and uniform horizontal flow is parallel to the axes of the convecting cells. One further paper which is of some interest in the present context is that of Rees [14], which considers the presence of viscous dissipation in a sidewall-heated square cavity. In that work, the author showed that the "velocity-squared" model for viscous dissipation in a porous medium, the form which is generally used and which was derived by Breugem and Ree [15] using a Representative Elementary Volume analysis, cannot be sustained by having values for the Gebhart number that are too large. In [14], it was found that, when Gebhart number increases beyond a critical value (which depends on the Rayleigh number but which is above unity), then there exists a region within the interior of the cavity that has temperatures that are larger than what is imposed on the boundary, a situation which is unphysical. Thus, large Gebhart numbers induce a physical situation that cannot occur in reality, and this concurs with the following statement by Celli et al. [12], namely, that: the condition Ge=1 is far beyond any conceivable application, except for some geophysical cases.
In the present paper, the Gebhart number is assumed to be asymptotically small, and its magnitude is used as the expansion parameter for the weakly nonlinear analysis.
It is shown that a hexagonal pattern of convection is preferred when convection first appears, and that this is due to the presence of viscous dissipation which introduces extra forcing terms at third order in the weakly nonlinear theory. At higher Rayleigh numbers within the weakly nonlinear regime, rolls are re-established as the favoured pattern, although there is an intermediate regime where both patterns are stable to small disturbances. It is found that, as the Forchheimer parameter increases from zero, hexagonal cells lose their characteristic subcriticality, but they increase their range of stability to small disturbances.
The analysis which is followed here is based heavily on the general framework for the weakly nonlinear theory of Darcy-Bénard convection, which is outlined in Rees [16], and the reader is referred to those lecture notes for further details.

Governing Equations
The equations governing Boussinesq convection in a porous medium subject to Darcy's law are given by ∂u ∂x where K is the permeability, g gravity, µ the fluid viscosity, k the thermal conductivity of the saturated medium and β the coefficient of cubical expansion. In addition, ρ is density and C p the specific heat, with the subscripts m and f referring, respectively, to the porous matrix and the saturating fluid. In the above, x and y are the horizontal coordinates, while z is the vertical coordinate with u, v and w being the respective Darcy velocities. In addition, T is the temperature of the saturated medium while T min is the reference temperature. Lengths are nondimensionalised with respect to d, the height of the layer, velocities with respect to k/d(ρC p ) f , time with respect to d 2 (ρC p ) m /k, and temperature using T = T min + (T max − T min )θ. Therefore, Equations (1)-(3) become ∂u ∂x + ∂v ∂y The viscous dissipation parameter, , is given by where are the Gebhart and Darcy-Rayleigh numbers, respectively. Physically, will be very small, and we will use this as the expansion parameter for the weakly nonlinear theory.
We may now eliminate the velocities from Equations (4) and (6) by using Equation (5). The basic pressure and temperature profiles corresponding to the steady conduction regime are We now perturb about this basic state by making the following substitutions: On dropping the tildes, the resulting nonlinear perturbation equations are and These equations are to be solved subject to the homogeneous boundary conditions, in a horizontally unbounded domain.

Weakly Nonlinear Analysis
Weakly nonlinear analysis is undertaken in the usual way by writing where the viscous dissipation parameter, , satisfies 0 < 1. The timescale must be slow, since we are close to the onset of instability, and therefore we set The factor of 1 2 in Equation (16) is an a posteriori numerical adjustment to maximise the number of unit coefficients in the final amplitude equations.
At O( ) in the expansion, the appropriate equations are given by These homogeneous equations admit roll solutions of any orientation, but at a critical value of R 0 , which depends on the wavenumber of the rolls. At this stage of the analysis, Equation (17) represents the linearised stability equations for the classical Darcy-Bénard problem, the neutral curve for which is given by R 0 = (a 2 + π 2 ) 2 /a 2 , where a is the wavenumber. The value of R 0 is minimised when a = π, and, therefore, we will concentrate on this value, for which R 0 = 4π 2 .
If we were to take as the solution of Equation (17) that of a single roll, then the nonlinear interaction caused by viscous dissipation terms in Equations (11) and (12) would not alter the final amplitude equation from that which arises when viscous dissipation is absent. However, the interaction of two rolls at 60 • to each other results in a third roll which is at 60 • to both of these. Thus, a quadratic nonlinearity occurs in the resulting set of three amplitude equations. Therefore, we take as the solution of Equation (17) the trio of rolls: where A, B and C are functions of the slow timescale τ, and the coefficients in Equations (18) and (19) have also been chosen to obtain as many unit coefficients as is possible in the amplitude equations obtained below.
At O( 2 ), the full equations are given by and the inhomogeneous terms expand to many lines and are therefore omitted for the sake of brevity. The solutions for p 2 and θ 2 are given by , the governing equations are given by Many of the inhomogeneous terms in Equations (24) and (25) are proportional to the O( ) eigensolutions given in Equations (18) and (19), and it is necessary to apply solvability conditions as given below. The final three terms on the right-hand side of Equation (25) correspond to the presence of viscous dissipation.
We take as our solvability condition where R 1 and R 2 are terms drawn from the right-hand sides of (24) and (25), respectively, and wherê p andθ are eigensolutions taken from one of the O( ) solutions in Equations (18) and (19). Then, the integration given in Equation (26) takes place over a region in which hexagonal patterns are periodic in both horizontal directions. A similar procedure takes place for the other two roll directions. The satisfaction of the integral constraint given in Equation (19) ensures that the inhomogeneous terms in Equations (24) and (25) do not have components that are proportional to the O( ) eigensolutions given in Equations (18) and (19). Therefore, Equations (24) and (25) may be solved if required; otherwise, Equations (18) and (19) do not possess solutions. Following this process for all three roll directions yields the following trio of amplitude equations for A, B and C: where Ω and α are given by, In general, the value of Ω depends on the relative orientation of the two rolls; Rees [16] shows that the general expression is when the relative orientation of two rolls is φ.

Analysis of the Amplitude Equations
We are interested in the formation of hexagonal cells, and therefore we obtain the following steady solution of Equations (27)-(29), which defines the value of E, the steady equilibrium solution. The plus sign corresponds to the upper branch while the minus sign corresponds to the lower branch; see Figure 1. This solution may be taken as being real, without loss of generality. The bifurcation diagrams corresponding to this solution and that of a single roll solution, e.g., A = R 1/2 2 , B = C = 0, are both shown in Figure 1. It is well-known that hexagonal solutions bifurcate from the zero solution via a transcritical bifurcation, and subcritical convection occurs at values of R 2 that are negative. In this case, hexagonal solutions exist when the following inequality is satisfied, as compared with rolls which exist only when R 2 > 0. We now describe very briefly a straightforward stability analysis of hexagonal solutions. If we perturb the solution given in Equation (32) about A = B = C = E using A = E + δA and so on, then we obtain the following linearised perturbation equations for δA, δB and δC: Given their symmetry, these equations may be simplified by setting δB = δC, on noting that it is not possible to set δB = −δC since then Equations (35) and (36) become inconsistent. Setting both δA and δB to be imaginary and proportional to exp λτ, where λ is real, results in two possible values of λ, namely, 0 and −3αE. Respectively, these mean that the whole of the solution curve is neutrally stable with respect to disturbances in phase, while that part of the lower branch which corresponds to negative values of E is unstable. Repetition of the same procedure with δA and δB taken to be real yields either In Equation (37), λ is positive only when positive values of E correspond to points below the turning point on the solution curve, while, in Equation (38), λ is positive when Therefore, we conclude that hexagons exist and are stable in the range and for points on the upper branch of the solution curve. Points on the lower branch are always unstable.

The Stability of Rolls
We turn now to determining the range of values of R 2 within which rolls are stable. Omitting the details, we may perturb about the single roll solution, A = R 1/2 2 , B = C = 0 in the same way as above. We obtain the perturbation equations, Equation (42) only admits decaying solutions when δA is real, and is again neutrally stable to disturbances in phase (i.e., when δA is imaginary). Therefore, Equations (43) and (44) must be considered to determine neutral stability. We find that the single roll solution is stable only when Summarising the above results, we find that hexagons exist when R 2 ≥ −151.621, and rolls exist when R 2 ≥ 0. Hexagons are stable when R 2 ≤ 26,633.8 and rolls are stable when R 2 ≥ 378.693. Thus, there is an intermediate regime within which both rolls and hexagons are stable to infinitesimal disturbances. Although we do not prove it here, there is a mixed mode branch of unstable solutions which bifurcates away from the roll branch at R 2 = 378.693 and then joins the branch of hexagonal solutions at R 2 = 26,633.8 as a subcritical bifurcation.

The Effect of Forchheimer Terms on Hexagonal Cells
We now consider the additional effect of form drag, as modelled by the Forchheimer terms. Following Rees [17], the non-dimensional equations that govern those cases where the quadratic Forchheimer terms are significant are given by ∂u ∂x where The constantK, which appears in Equation (51), is quoted by many papers to be given by Ergun's relation,K = 1.75L/150(1 − Φ), (Nield and Bejan [3]), where L is a lengthscale associated with the microstructure of the medium, and Φ is the effective porosity. Again, following Rees [17], we need to rescale G according to so that form drag effects are felt at O( 3 ), rather than at O( 2 ). The value, G * , is termed the Forchheimer parameter.
The expansion proceeds as before and with exactly the same solutions up to O( 2 ). At O( 3 ), the presence of form drag yields further inhomogeneous terms that must be included into the solvability condition. The O( 3 ) equations are now given by In Equations (53) and (54), the values u 1 , v 1 and w 1 are the velocities corresponding to the O( ) pressure distribution, p 1 , and the value q 1 is the O( ) velocity, The separate application of the solvability condition (26) for all three roll directions now yields the following amplitude equations: The final coefficient in these amplitude equations is given by The value of this integral depends on the relative magnitudes of A, B and C, and cannot be evaluated analytically. However, for steady rolls and hexagons, we have hexagons: The numerical value given in Equation (60) is consistent with that given in Rees [17]. If one were to wish to follow the time-evolution of A, B and C, then it would be necessary to evaluate F at every timestep.
We may now consider the effect of the Forchheimer terms on the amplitude on hexagonal convection. If we set A, B and C to all be equal to E, which is real, and assume that the flow is steady, then E satisfies the equation where we have set γ = 119.4882G * for simplicity of presentation. In this rather unusual amplitude equation, it is clear that E = 0 remains a possible solution, but the nonzero solutions now satisfy The solution curves for E fall into four different cases depending on the magnitudes of γ(>0) and α, and these are sketched in Figure 2. When γ = 0 form drag effects are absent, the flows considered earlier apply. When 0 < γ < α, there remains a subcritical flow on the upper branch of the solution curve, but the degree of subcriticality has now been reduced. In this range of values of γ, it is shown easily that the smallest value of R 2 on the upper solution branch is given by at which point the amplitude of convection is given by In this range of values of γ, the lower branch becomes weaker in the sense that |E| decreases when γ increases for any chosen value of R 2 ; this is a consequence of form drag reducing the efficacy of buoyancy forces to cause flow. There is also a discontinuous change of slope of the lower branch at At the point where γ = α, the upper branch becomes a parabola, which is characteristic of a standard supercritical pitchfork bifurcation. Therefore, subcritical flow does not exist when G * ≥ α/119.4882, or, equivalently, when Ge is small: As γ continues to increase above α, the amplitudes of convection at any chosen positive value of R 2 continues to decrease, and the onset of convection is via an asymmetric sharp-nosed bifurcation.
The lower branch remains unstable for small disturbances while the upper branch is stable initially as R 2 increases above zero. The equivalent situation for rolls, given in [13], is that the bifurcation is also sharp-nosed, but the bifurcation diagram remains symmetric.

Effect of the Forchheimer Terms on Stability
Finally, we present information on how the presence of form-drag affects the stability boundaries for hexagons and rolls. We omit all of the details and find that hexagons are unstable when while rolls are stable when Thus, hexagons increase and rolls decrease their respective ranges of stability when form drag becomes increasingly important.

Conclusions
We have sought to determine the combined effect of viscous dissipation and form drag on Darcy-Bénard convection. Viscous dissipation serves to cause the stable pattern of convection to change from rolls to hexagons for a range of values of the Darcy-Rayleigh number, which includes the critical value, Ra = 4π 2 , or, equivalently, R 2 = 0.
We have determined the range of stability of these patterns and compared them with that for rolls, which eventually re-emerge as the stable pattern as the Darcy-Rayleigh number increases. We find that form drag alters the shape of the solution curves to such an extent that, when it is sufficiently strong, there is no longer any subcritical convection, which is the ubiquitous characteristic of hexagonal patterns. Form drag also increases the size of the region within which hexagons are stable and that region within which rolls are unstable.
It is in the nature of weakly nonlinear theory that the phenomena it predicts take place close to the onset of convection. It provides an accurate indicator of the qualitative behaviour that will arise when governing parameters take larger values. In the present paper, we have considered both the viscous dissipation and form drag to be weak, and therefore the detailed behaviour of hexagonal cells which have been described take place asymptotically close to Ra = 4π 2 . Should the Gebhart number not be asymptotically small, then hexagonal convection will arise at values of the Darcy-Rayleigh number, which are at O(1) distances from the critical value, and therefore they will be seen easily in practice. In addition, if form drag effects are also strong, then the loss of subcriticality shown in Figure 2 will also be visible practically. The message from our analysis is that there will be O(1) values of both the Gebhart and Forchheimer parameters that also obey the qualitative behaviours which we have described here.
It is also worth stating that other agencies will alter the symmetry of the basic temperature field to cause hexagons to arise preferentially-this includes a temperature-dependent thermal diffusivity and weak vertical throughflow as first considered by Wooding [18] and Sutton [19]. However, the manner in which the weakly nonlinear theory progresses suggests that any changes to the present analysis will only be quantitative, and, in these cases, the general qualitative behaviour will be identical to that presented here.
Finally, we note that Rebhi et al. [20], who have studied the onset and nonlinear development of double-diffusive convection in a porous cavity, have identified regimes within which convection is subcritical in the sense that it is possible for nonlinear flow to exist at values of the Darcy-Rayleigh number that are below that given by linearised stability theory. The chief interest of this paper in the present context is the observation that an increasingly strong form-drag also removes the subcriticality of the cavity flows. It is difficult to know whether this qualitative similarity between two very different flows has a great significance, for the mechanisms which cause subcritical convection are quite different in the respective cases.