## Abstract

A model is presented for thermal convection in an inclined layer of porous material when the medium has a bidispersive structure. Thus, there are the usual macropores which are full of a fluid, but there are also a system of micropores full of the same fluid. The model we employ is a modification of the one proposed by Nield & Kuznetsov (2006 *Int. J. Heat Mass Transf.* **49**, 3068–3074. (doi:10.1016/j.ijheatmasstransfer.2006.02.008)), although we consider a single temperature field only.

## 1. Introduction

Thermal convection in an inclined layer is a subject that has attracted much attention. For the situation where the layer is full of a Newtonian fluid a comprehensive analysis is provided by Chen & Pearlstein [1], who also give a lucid historical review of work prior to their article. Recently, many authors have been investigating thermal convection in a fluid-saturated, inclined porous layer, or even a vertical layer, see e.g. Barletta [2], Barletta & Celli [3], Barletta & Storesletten [4] and Barletta & Rees [5]; the last-mentioned contribution contains a good review of earlier work.

In this article, we shall commence an analysis of thermal convection in a fluid-saturated porous layer but when the porous material is of bidispersive type. Bidispersive porous materials were introduced theoretically in the late 1990s, e.g. Nield [6,7], Nield & Bejan [8], Nield & Kuznetsov [9] and chapter 13 of Straughan [10]. See also the recent paper by Nield [11]. A bidispersive porous medium is one where there are pores of the usual type which are referred to as macropores. Associated with the macropores is the porosity *ϕ* which is the ratio of the volume of the macropores to the total volume of saturated porous material.

However, in the solid skeleton, there are much smaller micropores, e.g. the picture given in Nield & Kuznetsov [9], p. 3069. Attached to the micropores is a porosity *ϵ* which is the ratio of the volume occupied by the micropores to the volume of the porous body that remains once the macropores are removed. Thus, the fraction of volume occupied by the micropores is *ϵ*(1−*ϕ*), whereas the fraction occupied by the solid skeleton is (1−*ϵ*)(1−*ϕ*). Nield & Kuznetsov [9] gave a comprehensive analysis of linear instability for thermal convection in a bidispersive porous material when the layer is horizontal. They allowed for different velocities *T*^{f} and *T*^{p} in the macro and micropores. Our goal is here to assess the effect of induced velocity owing to inclination and gravity upon instability, and so we here restrict attention to a single temperature *T*, although we do follow Nield & Kuznetsov [9] and allow different macro- and micropore velocities

We should point out that bidispersive porous materials (or double porosity materials) are very important and are the subject of intense recent research owing to the many real-life application areas. For example, bidispersive porous media are crucial to understanding the physics of underground drinking water supplies [12,13]; they are also crucial to understanding the highly controversial topic of hydraulic fracturing (or ‘fracking’) in which rocks underground are deliberately vibrated to release trapped natural gas [14,15]. Other application areas are given in Straughan [10]. A further area of particular interest for thermal convection in an inclined plane is to landslides [16] or more generally land deformation involving thermal gradients [17,18].

In §2, we introduce the basic equations for thermal convection in a bidispersive porous medium. In particular, we derive the energy balance law when a single temperature is employed. In the section following, we derive the basic solution for thermal convection in an inclined layer. The section following that analyses linear instability in detail. The penultimate section presents some global stability results when the fully nonlinear theory is employed. The final section is reserved for conclusions.

## 2. Basic equations

Let *T*. Thus, setting *T*^{f}=*T*^{p}=*T* in equations (6)–(10) of Nield & Kuznetsov [9] and omitting the Brinkman terms, the equations for

In these equations, *μ* is the fluid viscosity, *K* refers to permeability, *ζ* is an interaction coefficient, *p* refers to the pressures and sub/superscript ‘f’ and ‘p’ refer to macro and micropore effects. In addition, *T*_{F} is a reference temperature where the buoyancy term involves a linear density of form
*g*_{i} to denote the gravity vector.

We suppose that the porous material occupies an inclined plane with angle from the vertical *δ*, where 0≤*δ*≤*π*/2 radians, cf. the geometrical configuration of Chen & Pearlstein [1]. Thus, in our case,
*x* is measured in the longitudinal direction along the layer, *z* is orthogonal, with the boundaries being at *z*=0 and *z*=*d*, and **g** is gravity.

To derive the energy balance equation governing the behaviour of the temperature field, we generalize the procedure of Joseph [19], see also Straughan [20], p. 47 and [10], p. 16 . Let *V* be a representative elementary volume, REV, and denote by *s*, *f* and *p* the properties of the solid skeleton, macropores and micropores, respectively. Let *V* ^{f}_{i} and *V* ^{p}_{i} be the actual velocities in the macro- and micropores. Then, we may write equations for the separate temperature fields *T*^{s}, *T*^{f} and *T*^{p} as
*ρ* and *c* refer to density and specific heat at constant pressure, *κ* is the thermal conductivity and *h* is a temperature interaction coefficient. Throughout, we employ standard indicial notation, and Δ is the *Laplace* operator. We have followed Nield & Kuznetsov [9] and employed interactions between *T*^{f} and *T*^{p}. We could also include interactions between *T*^{s} and *T*^{f} or *T*^{s} and *T*^{p}.

In the next step, we multiply (2.2)_{1} by (1−*ϕ*)(1−*ϵ*), (2.2)_{2} by *ϕ* and (2.2)_{3} by *ϵ*(1−*ϕ*), and add the results, setting *T*^{s}=*T*^{f}=*T*^{p}=*T*.

In the REV, we thus have, recalling

We have followed Nield & Kuznetsov [9] and allowed (*ρc*)_{f} and (*ρc*)_{p} to be different. Because Nield & Kuznetsov dealt with different macro- and micropore temperatures, it is possible (*ρc*)_{f} and (*ρc*)_{p} will be different. However, in the present case of single temperature, we expect that as the macro- and micropores are composed of the same fluid (*ρc*)_{f} and (*ρc*)_{p} will be the same. Thus, in our case, we may simplify equations (2.1) and (2.2), and we work with the system
*T* in a Darcy bidispersive porous medium with a single temperature may be taken to be (2.3).

## 3. Perturbation equations

From equation (2.3), we seek a basic solution of the form

for constants *T*_{L} and *T*_{U}. We consider the two cases of heating from above, *T*_{U}>*T*_{L} and heating from below, *T*_{L}>*T*_{U}. In addition, we normalize the solution by requiring
*z*≤*d* into *z*∈[0,1], we find the basic solution in whose stability we are interested is
*T*|=|*T*_{L}−*T*_{U}|.

The basic velocities have form
*Ξ* is given by

We now introduce perturbations and write

We non-dimensionalize with length scale *d*, time scale
*K*_{r}=*K*_{p}/*K*_{f}, Rayleigh number
*σ*_{f}=*ζK*_{f}/*μ* and *N*_{k},*K*_{N} are non-dimensional parameters of form
*K*_{N} is a sort of thermal Reynolds number.

## 4. Linear instability

To proceed, we linearize equations (3.1), then remove the pressures *π*^{f} and *π*^{p} by taking *curlcurl* of equations (3.1)_{1,2}. Then, in terms of *w*^{f}, *w*^{p} and *θ*, we obtain equations and then introduce normal modes by writing *f*=*f*(*z*) e^{i(ax+by)+ct}, where *f* is *w*^{f}, *w*^{p} and *θ*, to obtain
*D*=*d*/*dz*.

Equations (4.1)–(4.2) are solved numerically by a Chebyshev tau method [21] with boundary conditions

## 5. Nonlinear stability

In this section, we include a global nonlinear stability result when the inclined layer is heated from above.

We commence with equations (3.1) with the top signs, so that we have

We suppose that
*θ*, *π*^{f} and *π*^{p} satisfy a plane tiling periodicity in *x* and *y*, such that the solution has a periodicity cell *V* .

Multiply (5.1)_{1} by _{2} by *V* . Likewise multiply (5.1)_{3} by *V* . After some integrations by parts and use of the boundary conditions, we add the results to obtain the equation

In (5.2), *L*^{2}(*V*). It is worth observing that the term involving *θ*_{,x} integrates to zero as do the cubic nonlinear terms which arise.

From (5.2), we obtain
_{1} the first eigenvalue in the membrane problem for *V* (in fact λ_{1}=*π*^{2}) and then we use the arithmetic–geometric mean inequality to find
*ϵ*_{1}>0,*ϵ*_{2}>0 are at our disposal.

We now select *θ*∥ decays at least exponentially in time. Thus, (5.5) represents a global nonlinear stability threshold.

To establish decay of **u**^{f} and **u**^{p} under (5.5), we employ the relations obtained above to find after using arithmetic–geometric mean inequality
*θ*∥ decays exponentially, (5.6) demonstrates ∥**u**^{f}∥ and ∥**u**^{p}∥ likewise have at least exponentially decay in *t*.

It is worth observing that as *δ*→*π*/2 the right-hand side of inequality (5.5) tends to

We observe that in the case of an inclined layer heated from below, it is possible to obtain a global nonlinear stability result following the estimates given in [22]. A global stability condition that may be obtained is the following

## 6. Numerical results and conclusions

In this section, we present the results for the numerical solution of (4.1)–(4.3), when the layer is heated from below. We employ the Chebyshev tau method and allow for minimization in both *a* and *b* to find the minimum value *Ra*=*R*^{2}. In this section, we take *σ*_{f}=1. We have employed other values of *σ*_{f} and the variation in the results is not significant. However, the bidispersive parameter *K*_{r} plays a strong role as we report below.

Figures 1 and 2 show the critical surfaces in the (*a*,*b*,*Ra*) space, i.e. the surface defined implicitly by *c*(*a*,*b*,*Ra*)=0 (*c* is the time evolution exponent, see equation (4.2)), for various angles of inclination (*φ*=*π*/2−*δ*). These surfaces are analogous to the critical curves, such as those of figure 4, when we have to consider the two wavenumbers *a*,*b* instead of just *a*. We see that for all angles of inclination, once *φ*>0, the lowest value for *Ra* is found when *a*=0, *b*=*π* (the minimum of *Ra* on the black line on the plane *a*=0), so the most unstable convection cells are longitudinal rolls. This is what we might expect, because the same is true for inclined convection in the classical problem employing Darcy’s law, cf. e.g. Rees *et al.* [23], p. 16. In figure 1, the surface of instability touches the plane *b*=0. In figure 2 with *φ*=17°, there is still an intersection with *b*=0, but the intersection is reduced. Once *φ*=20°), there is no intersection at all with the plane *b*=0, so all perturbations in the form of transverse rolls are stable.

Figure 3 shows the critical Rayleigh number values versus the inclination angle *φ* for various values of *K*_{r}=*K*_{p}/*K*_{f}. The transverse mode curves finish abruptly as shown, because there is a critical angle where any instability owing to transverse rolls disappears. This is in complete agreement with the work of Rees & Bassom [24] and the classical inclined Darcy convection problem. Figure 3 demonstrates that the critical Rayleigh number corresponding to longitudinal rolls coincides with that for transverse rolls when *φ*=0 but is lower immediately when *φ*>0 (*this means instability occurs via longitudinal rolls*). Figures 4 and 5 show the behaviour of *Ra* versus *a* for transverse rolls when *K*_{r}=10^{−4} for various angles *φ*. This behaviour is in agreement with that reported for the single porosity case by Rees & Bassom [24]. As is seen from figure 4, a disconnected and closed instability curve forms once *φ* exceeds 16.255°. Figure 5 shows how this curve behaves as *φ* increases beyond 17°. A numerical analysis shows that the upper unbounded part of the critical curve shown in figure 5 disappears for some *φ*<17.5°, whereas the lower closed curve persists at least for *φ*≈17.52°. The same qualitative behaviour is found for other values of *K*_{r}, so all the transverse mode curves in figure 3 are continuous until their end. We find (for *K*_{r}=10^{−4}) no instability by a transverse mode once

While longitudinal modes do yield the linear instability boundary, it may be that this problem resembles more problems of parallel shear flow, where more than the leading eigenvalue influences the global stability scenario, cf. Straughan [25], ch. 8.

To make a direct comparison with the work of Rees & Bassom [24], we take *K*_{r}=1 and then form the sum of equations (3.1)_{1} and (3.1)_{2}. This yields a system of equations of the form
*v*_{3}=*w*^{f}+*w*^{p}. Make the transformation *S*^{2}=*Ra*_{1}=2*R*^{2}=2 *Ra*. Thus, the Rayleigh number reported here is one half of the number reported by Rees & Bassom [24]. We find the angle when instability disappears when *K*_{r}=1 is ≈31.25°, which agrees with [24]. We do find that the bidispersive case is very different from the simple porosity case as figure 3 indicates. In fact, for longitudinal rolls equations (4.1)–(4.3) may be solved analytically to yield
*Z* the quantity (2+*K*_{r})/(1+5*K*_{r}). Then, the critical eigenfunctions *w*^{f}, *w*^{p} and *θ* are given by
*K*_{r}, is playing a key role in determining the forms of *w*^{f}, *w*^{p} and *θ* at the onset of instability. It is worth pointing out that we found no overstability numerically. By checking the error coefficient *τ*, we have found that any complex eigenvalues which arise are associated with spurious eigenvalues.

## Ethics

This work does not pose ethical issues.

## Data accessibility

This work has no supplementary data.

## Authors' contributions

B.S. proposed the model and drafted the introduction and the concluding section; P.F. carried out the linear instability computations and helped draft the concluding section; G.M. contributed to the nonlinear stability section. All authors gave final approval for publication.

## Competing interests

We have no competing interests.

## Funding

Research partially supported by: the University of Catania under an FIR contract; ‘Gruppo Nazionale della Fisica Matematica’ of the ‘Istituto Nazionale di Alta Matematica’

## Acknowledgements

We acknowledge the contribution of University of Catania and the INDAM.

- Received June 15, 2016.
- Accepted July 1, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.