1 Introduction
Fluid matter exhibits a remarkable tendency to build up patterns and organized dynamical structures. Different stages can appear in a fluid system as it moves further away from equilibrium (Gollup Reference Gollup1995; Kadanoff Reference Kadanoff2001): laminar convection ends up developing turbulence (Batchelor Reference Batchelor1982) that produces different characteristic patterns (Hunt & Durbin Reference Hunt and Durbin1999), such as plumes, swirls, eddies, vortices (the vortex inside Saturn’s hexagon is a beautiful example of an atmospheric stable vortex, see for instance the work by Godfrey (Reference Godfrey1990)) and, eventually, spatio-temporal chaos (Egolf et al. Reference Egolf, Melnikov, Pesch and Ecke2000). Furthermore, extensive spatio-temporal chaos may render equilibrium-like states back again at larger time/length scales (Egolf Reference Egolf2000).
Since the seminal works by Bénard (Reference Bénard1900) and Lord Rayleigh (1916), the fluid flow in closed circuits (convection), caused by the presence of temperature inhomogeneities, has likely been one of the most studied problems in science (the works by Cross & Hohenberg (Reference Cross and Hohenberg1993), Bodenschatz, Pesch & Ahlers (Reference Bodenschatz, Pesch and Ahlers2000) and Mutabazi, Wesfreid & Guyon (Reference Mutabazi, Wesfreid and Guyon2006) are good reviews on the subject). The phenomenon is well known to be ubiquitous in nature, including biological systems. However, let us describe it again, at its simplest: we consider a real, experimental system where there is a gravity force keeping a fluid layer at rest at temperature
$T$
and with two horizontal and limiting surfaces. The upper one (in the sense of gravity) is either free or in contact with a solid surface. Then, by means of some kind of temperature source at higher temperature
$T_{0}>T$
, the fluid is heated from below. The difference
$T_{0}-T$
is gradually increased, but the fluid remains static. However, when a critical temperature gradient is reached, fluid motion is set up, shaping regular patterns in all of the fluid volume (Bénard Reference Bénard1900). According to classical theory, the Rayleigh number is the convection control parameter in this kind of problem. This dimensionless parameter is usually defined as
$Ra\equiv \unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}Tgh^{3}/\unicode[STIX]{x1D705}\unicode[STIX]{x1D708}$
, where
$\unicode[STIX]{x1D6FC}\equiv (1/V)(\unicode[STIX]{x2202}T/\unicode[STIX]{x2202}V)_{p}$
is the fluid expansion (
$V$
is the fluid volume and
$T$
and
$p$
are the temperature and hydrostatic pressure, respectively) coefficient,
$\unicode[STIX]{x0394}T$
is the boundaries temperature difference,
$g$
is gravity acceleration,
$h$
is the system width and
$\unicode[STIX]{x1D705}$
and
$\unicode[STIX]{x1D708}$
are the thermal conductivity and kinematic viscosity transport coefficients, respectively. In fact, linear theory predicts a critical value
$Ra_{c}=657.5$
and
$Ra_{c}=1708$
for the free surface and the closed-on-top system cases, respectively (Mutabazi et al.
Reference Mutabazi, Wesfreid and Guyon2006). The theoretical treatment by Lord Rayleigh (1916) relied on the work by Boussinesq (Reference Boussinesq1903) (see also the book edited by Mutabazi et al. (Reference Mutabazi, Wesfreid and Guyon2006) for a more recent review), who defined the relevant contributions for this convection in the fluid balance equations. These balance equations, as is known, were worked out by Navier and Stokes only a few years before (Batchelor Reference Batchelor1967).
Let us recall that a more generic concept of fluid also involves systems where the particles are not necessarily microscopic, i.e. the particles can be macroscopic, when their typical size is greater than
$1~\unicode[STIX]{x03BC}\text{m}$
(Bagnold Reference Bagnold1954). In fact, the dynamical properties of a set of rigid macroscopic particles in a high state of agitation was elucidated as a subject of the theory of fluids a long time ago by Reynolds (Reference Reynolds1885). However, for macroscopic particles, the kinetic energy will be partially transferred, upon collision, to the lower (smaller length scales) dynamics levels, never coming back to the upper granular level. For instance, it may be transferred into thermal movement of the molecules that are the constituents of the disk material (Andreotti, Forterre & Pouliquen Reference Andreotti, Forterre and Pouliquen2013). Thus, unless the system gets an energy input from some kind of source, it will evolve by continuously decreasing its total kinetic energy, i.e. lowering the system granular temperature (Kanatani Reference Kanatani1979). This temperature decay rate was calculated, for a homogeneous and low-density granular system (i.e. a homogeneous granular gas), by Haff (Reference Haff1983).
Nevertheless, when excited by some persistent external action, stable granular gas systems are found spontaneously in nature, for instance, in sand storms (Bagnold Reference Bagnold1954), and also in laboratory experiments, where air flow (Losert et al. Reference Losert, Bocquet, Lubensky and Gollub2000) or mechanical vibration (Olafsen & Urbach Reference Olafsen and Urbach1998; Puglisi et al. Reference Puglisi, Gnoli, Gradenigo, Sarracino and Villamaina2012) may be used as energy inputs. Under these conditions, the granular gas can develop steady laminar flows (Vega Reyes & Urbach Reference Vega Reyes and Urbach2009). Unfortunately, the hydrodynamics of granular fluids is not in the same stage of development as it is for molecular fluids (Puglisi Reference Puglisi2015). For instance, the corresponding hydrodynamic theory for thermal convection in a granular gas was developed only very recently (see, for instance, the work by Khain & Meerson (Reference Khain and Meerson2003)) and only in the case of the academic problem of horizontal (or inclined) infinite walls. However, real systems are obviously finite. In the case of a gas heated by two horizontal walls, the simplest finite configuration considered is when the system is closed by adding vertical lateral walls. As is known, finite-size effects have an impact on both the critical Rayleigh number and the convection scenario in molecular fluid changes if a cold vertical wall is present (Daniels Reference Daniels1977; Hall & Walton Reference Hall and Walton1977; Mutabazi et al. Reference Mutabazi, Wesfreid and Guyon2006).
The lack of a theoretical analysis on the effects of finite-size systems in granular convection theory (see, for instance, the works by Forterre & Pouliquen (Reference Forterre and Pouliquen2003), Khain & Meerson (Reference Khain and Meerson2003), Eshuis et al. (Reference Eshuis, van der Weele, van der Meer and Lohse2007), Bougie (Reference Bougie2010), Eshuis et al. (Reference Eshuis, van der Meer, Alam, van Gerner, van der Weele and Lohse2010) and others) may have hindered and/or rendered impossible a complete interpretation of a part of the previous results on granular dynamics laboratory experiments. Furthermore, as we will see in some experimental works, the observed granular convection (Wildman, Huntley & Parker Reference Wildman, Huntley and Parker2001a ,Reference Wildman, Huntley and Parker b ; Risso et al. Reference Risso, Soto, Godoy and Cordero2005; Eshuis et al. Reference Eshuis, van der Meer, Alam, van Gerner, van der Weele and Lohse2010; Windows-Yule, Rivas & Parker Reference Windows-Yule, Rivas and Parker2013; Pontuale et al. Reference Pontuale, Gnoli, Vega Reyes and Puglisi2016) should be either exclusively or partly due to sidewall energy sink, and not of the Rayleigh–Bénard type. Of course, although the no-sidewalls theoretical approach may be accurate when bulk convection is present (Khain & Meerson Reference Khain and Meerson2003), in the cases where the convection is caused only by the sidewalls energy sink, we may expect the convection properties to be very different. Note that in an enclosed granular gas, the lateral energy sink should always be present, since wall–particle collisions are inherently inelastic. This implies that the present analysis should be relevant for many granular convection experiments. Furthermore, as we will see, lateral-wall effects are also more substantial for the granular gas than for the molecular fluid.
2 Description of the system and the problem
Let us consider a system consisting of a large set of circular particles (disks) in a two-dimensional (2D) system. The particles are identical inelastic smooth hard disks with mass
$m$
and diameter
$\unicode[STIX]{x1D70E}$
. The hard-collision model works reasonably well at an experimental level for a variety of materials (Foerster et al.
Reference Foerster, Louge, Chang and Allis1994). We use here this model in the smooth particle approximation, i.e. we neglect the effects of sliding and friction in the collision. Under the smooth hard particle model for collisions, the fraction of kinetic energy loss after collision is characterized by a constant parameter called the coefficient of normal restitution
$\unicode[STIX]{x1D6FC}$
, not to be confused with the expansion coefficient, usually also denoted by
$\unicode[STIX]{x1D6FC}$
:
$$\begin{eqnarray}\displaystyle \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{v}_{12}^{\prime }=-\unicode[STIX]{x1D6FC}\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{v}_{12}, & & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{v}_{12},\boldsymbol{v}_{12}^{\prime }$
are the collision-pair contact velocities before and after collision, respectively (Foerster et al.
Reference Foerster, Louge, Chang and Allis1994).
The system is under the action of a constant gravitational field
$\boldsymbol{g}=-g\hat{\boldsymbol{e}}_{y}$
. We will also assume that the system has low particle density (
$n$
) everywhere at all times. Therefore, our fluid is a granular gas. Collisions are instantaneous (in the sense that the contact time is very short compared with the average time between collisions, Chapman & Cowling Reference Chapman and Cowling1970) and occur only between two particles. Since particles are inelastic, our theory should take into account the inelastic cooling term.
The system is provided with either two (top and bottom) or just one (bottom, in the sense of gravity) horizontal walls (i.e. perpendicular to gravity), these being provided with energy sources. In addition, our granular gas is caged in a finite rectangular region by two inert vertical walls (we call them lateral walls or sidewalls). Sidewalls–particle collisions are inherently inelastic, the degree of inelasticity of these collisions being characterized by a coefficient of normal restitution
$\unicode[STIX]{x1D6FC}_{w}$
that is, in general, different from the one for particle–particle collisions,
$\unicode[STIX]{x1D6FC}$
. See figure 1 for a graphical description.
Simple sketch of the system. The system is heated from the horizontal walls (at
$y=\pm h/2$
). If the lateral walls (at
$x=\pm L/2$
) act as energy surface sinks, bidimensional flow occurs.
We denote the single-particle velocity distribution function as
$f(\boldsymbol{r},\boldsymbol{v}|t)$
with
$\boldsymbol{r},\boldsymbol{v}$
being the particle position and velocity, respectively, functions of time (
$t$
). The first three velocity moments of the distribution function,
$n(\boldsymbol{r},t)=\int \text{d}\boldsymbol{v}f(\boldsymbol{r},\boldsymbol{v}|t)$
,
$\boldsymbol{u}(\boldsymbol{r},t)=(1/n)\int \text{d}\boldsymbol{v}f(\boldsymbol{r},\boldsymbol{v}|t)\boldsymbol{v}$
and
$T(\boldsymbol{r},t)=(1/dn)\int \text{d}\boldsymbol{v}f(\boldsymbol{r},\boldsymbol{v}|t)mV^{2}$
, define the average fields’ particle density (
$n$
), flow velocity (
$\boldsymbol{u}$
) and temperature (
$T$
), respectively. Here,
$\boldsymbol{V}=\boldsymbol{v}-\boldsymbol{u}$
and
$d$
is the system dimension. In this work, we consider only
$d=2$
(and that is why the particles are necessarily flat).
For a granular gas, molecular chaos (i.e. particle velocities are not statistically correlated) also occurs in most practical situations (Prevost, Egolf & Urbach Reference Prevost, Egolf and Urbach2002; Baxter & Olafsen Reference Baxter and Olafsen2007). Therefore, the kinetic Boltzmann equation (Chapman & Cowling Reference Chapman and Cowling1970), may also be used to describe granular gases (Brey et al. Reference Brey, Dufty, Kim and Santos1998; Dufty Reference Dufty2001). The general balance equations that follow from the inelastic Boltzmann equation have the same form as that for molecular fluid except for the additional inelastic cooling term arising in the energy equation. They have the following form (Brey et al. Reference Brey, Dufty, Kim and Santos1998; Sela & Goldhirsch Reference Sela and Goldhirsch1998)
$$\begin{eqnarray}\displaystyle D_{t}n=-n\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u},\quad D_{t}\boldsymbol{u}=-\frac{1}{mn}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}+\boldsymbol{g},\quad D_{t}T+\unicode[STIX]{x1D701}T=-\frac{2}{dn}(\unicode[STIX]{x1D64B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{q}). & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
In the above equations,
$D_{t}\equiv \unicode[STIX]{x2202}_{t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$
is the material derivative (Batchelor Reference Batchelor1967), and
$\unicode[STIX]{x1D64B},\boldsymbol{q}$
are the moment and energy fluxes (stress tensor and heat flux), defined by
$\unicode[STIX]{x1D64B}=m\int \text{d}\boldsymbol{v}\boldsymbol{V}\boldsymbol{V}f(\boldsymbol{v})$
and
$\boldsymbol{q}=(m/2)\int \text{d}\boldsymbol{v}V^{2}\boldsymbol{V}f(\boldsymbol{v})$
, respectively.
As stated previously, note the new term
$\unicode[STIX]{x1D701}T$
in the energy equation of (2.2), where
$\unicode[STIX]{x1D701}$
represents the rate of kinetic energy loss, and is usually called the inelastic cooling rate. Accurate expressions of the cooling rate and the inelastic Boltzmann equation are well known and may be found elsewhere (we use here that worked out by Brey et al. (Reference Brey, Dufty, Kim and Santos1998)).
Let us note also that the set of balanced equations (2.2) is exact and always valid. However, to close the system of equations, we need to express the fluxes
$\unicode[STIX]{x1D64B},\boldsymbol{q}$
and the cooling rate
$\unicode[STIX]{x1D701}$
as functions of the average fields
$n,\boldsymbol{u},T$
. The hydrostatic pressure field
$p$
is defined by the equation of state for an ideal gas:
$p=nT$
. Starting out of the kinetic equation for the gas, this can only be done if the distribution function spatio-temporal dependence can be expressed through a functional dependence on the average fields, i.e. if the gas is in a normal state (Hilbert Reference Hilbert1912), where this is only true if the spatial gradients vary over distances greater than the mean free path (that is, the characteristic microscopic scale). In this case, it is usually said that there is scale separation (Goldhirsch Reference Goldhirsch2003). Henceforth, we assume this scale separation occurs at all situations considered for this work (see Vega Reyes & Urbach (Reference Vega Reyes and Urbach2009) for more details about the conditions for accuracy of this assumption in steady granular gas flows).
The boundary conditions come from the usual forms for temperature sources, no-slip velocity (Vega Reyes, Santos & Garzó Reference Vega Reyes, Santos and Garzó2010) and controlled pressure at the boundaries:
$T(y=+h/2)=T_{+}$
, (substituted by
$[\unicode[STIX]{x2202}T/\unicode[STIX]{x2202}y=0]_{y=+h/2}$
, with
$h\gg L$
, in the case of an open-on-top system),
$T(y=-h/2)=T_{-}$
,
$\boldsymbol{u}(x=\pm L/2)=\boldsymbol{u}(y=\pm h/2)=\mathbf{0}$
,
$p(y=-h/2)=p_{0}$
. For an enclosed granular gas, we also necessarily need to consider the dissipation at the lateral walls, as we explained previously. A condition for the horizontal derivative of temperature would suffice to account for an energy sink at the side walls (Hall & Walton Reference Hall and Walton1977). However, taking into account that the energy sink comes from wall–particle collision inelasticity, then it is more appropriate to assume a horizontal heat flux that is proportional to lateral wall–disks collisions’ degree of inelasticity,
$\propto (1-\unicode[STIX]{x1D6FC}_{w}^{2})$
(Johnson & Jackson Reference Johnson and Jackson1987),
$$\begin{eqnarray}\displaystyle q_{x}(x=\pm L/2)={\mathcal{A}}(\unicode[STIX]{x1D6FC}_{w})[pT^{3/2}]_{x\,=\,\pm L/2}, & & \displaystyle\end{eqnarray}$$
where
${\mathcal{A}}(\unicode[STIX]{x1D6FC}_{w})=(\unicode[STIX]{x03C0}/2)m(1-\unicode[STIX]{x1D6FC}_{w}^{2})$
is given by the dilute limit of the corresponding expression in the work by Nott et al. (Reference Nott, Alam, Agrawal, Jackson and Sundaresan1999). (In addition, condition
$p(x=-L/2)=p(x=+L/2)$
is also used implicitly, since in this work, we only consider identical lateral walls at both sides.) The detailed balance of fluxes across the boundaries is beyond the scope of this work, but for a more detailed analysis on realistic boundary conditions, the reader may refer to the work by Nott et al. (Reference Nott, Alam, Agrawal, Jackson and Sundaresan1999).
2.1 Navier–Stokes equations and transport coefficients
We assume that the spatial gradients are sufficiently small, which is true for steady laminar flows near the elastic limit (Vega Reyes & Urbach Reference Vega Reyes and Urbach2009). Therefore, we use the Navier–Stokes constitutive relations for the fluxes (Brey et al. Reference Brey, Dufty, Kim and Santos1998; Brey & Cubero Reference Brey, Cubero, Pöschel and Luding2001)
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D64B}=p\unicode[STIX]{x1D644}-\unicode[STIX]{x1D702}\left[\unicode[STIX]{x1D735}\boldsymbol{u}+\unicode[STIX]{x1D735}\boldsymbol{u}^{\dagger }-\frac{2}{d}\unicode[STIX]{x1D644}(\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u})\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{q}=-\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}T-\unicode[STIX]{x1D707}\unicode[STIX]{x1D735}n. & \displaystyle\end{eqnarray}$$
In (2.5) we can find the transport coefficients:
$\unicode[STIX]{x1D702}$
(viscosity),
$\unicode[STIX]{x1D705}$
(thermal conductivity), and
$\unicode[STIX]{x1D707}$
(thermal diffusivity). Note that
$\unicode[STIX]{x1D707}$
is a new coefficient that arises from inelastic particle collisions (Brey et al.
Reference Brey, Dufty, Kim and Santos1998). The transport coefficients for the granular gas have been calculated by several authors, with some variations in the theoretical approach. For instance, Sela & Goldhirsch (Reference Sela and Goldhirsch1998) performed a Chapman–Enskog-like power series expansion in terms of both the spatial gradients and inelasticity up to the Burnett order, but limited to quasi-elastic particles, whereas Brey et al. (Reference Brey, Dufty, Kim and Santos1998) perform the expansion only in the spatial gradients (and the theory is formally valid for all values of inelasticity). Previous works, such as the work by Jenkins & Savage (Reference Jenkins and Savage1983) and by Lun et al. (Reference Lun, Savage, Jeffrey and Chepurniy1984) obtain the granular gas transport coefficients only in the quasi-elastic limit. These theories will actually yield indistinguishable values of the Navier–Stokes transport coefficients for nearly elastic particles (the case of our interest in the present work).
The essential point to our problem is the scaling with temperature
$T$
and particle density
$n$
. This scaling for hard particles is (Chapman & Cowling Reference Chapman and Cowling1970; Brey et al.
Reference Brey, Dufty, Kim and Santos1998)
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}=\unicode[STIX]{x1D702}_{0}^{\ast }T^{1/2},\quad \unicode[STIX]{x1D705}=\unicode[STIX]{x1D705}_{0}^{\ast }T^{1/2},\quad \unicode[STIX]{x1D707}=\unicode[STIX]{x1D707}_{0}^{\ast }\frac{T^{3/2}}{n},\quad \unicode[STIX]{x1D701}=\unicode[STIX]{x1D701}_{0}^{\ast }\frac{p}{T^{1/2}}. & & \displaystyle\end{eqnarray}$$
The values of the coefficients for disks are
$\unicode[STIX]{x1D702}_{0}^{\ast }\equiv \unicode[STIX]{x1D702}^{\ast }(\unicode[STIX]{x1D6FC})\sqrt{m}/(2\unicode[STIX]{x1D70E}\sqrt{\unicode[STIX]{x03C0}}),\unicode[STIX]{x1D705}_{0}^{\ast }\equiv 2\unicode[STIX]{x1D705}^{\ast }(\unicode[STIX]{x1D6FC})/(\sqrt{\unicode[STIX]{x03C0}m\unicode[STIX]{x1D70E}}),\unicode[STIX]{x1D707}_{0}^{\ast }\equiv 2\unicode[STIX]{x1D707}^{\ast }(\unicode[STIX]{x1D6FC})/(\sqrt{\unicode[STIX]{x03C0}m\unicode[STIX]{x1D70E}}),\unicode[STIX]{x1D701}_{0}^{\ast }\equiv (2\unicode[STIX]{x1D70E}\sqrt{\unicode[STIX]{x03C0}/m})\unicode[STIX]{x1D701}^{\ast }(\unicode[STIX]{x1D6FC})$
the expressions’ dimensionless functions that depend on the coefficient of restitution can be found in the work by Brey & Cubero (Reference Brey, Cubero, Pöschel and Luding2001). We write them here for completeness
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}^{\ast }(\unicode[STIX]{x1D6FC})=\left[\unicode[STIX]{x1D708}_{1}^{\ast }(\unicode[STIX]{x1D6FC})-\frac{\unicode[STIX]{x1D701}^{\ast }(\unicode[STIX]{x1D6FC})}{2}\right]^{-1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D705}^{\ast }(\unicode[STIX]{x1D6FC})=\left[\unicode[STIX]{x1D708}_{2}^{\ast }(\unicode[STIX]{x1D6FC})-\frac{2d}{d-1}\unicode[STIX]{x1D701}^{\ast }(\unicode[STIX]{x1D6FC})\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}^{\ast }(\unicode[STIX]{x1D6FC})=2\unicode[STIX]{x1D701}^{\ast }(\unicode[STIX]{x1D6FC})\left[\unicode[STIX]{x1D705}^{\ast }\left(\unicode[STIX]{x1D6FC}+\frac{(d-1)c^{\ast }(\unicode[STIX]{x1D6FC})}{2d\unicode[STIX]{x1D701}^{\ast }(\unicode[STIX]{x1D6FC})}\right)\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D701}(\unicode[STIX]{x1D6FC}^{\ast })=\frac{2+d}{4d}(1-\unicode[STIX]{x1D6FC}^{2})\left[1+\frac{3}{32}c^{\ast }(\unicode[STIX]{x1D6FC})\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D708}_{1}^{\ast }(\unicode[STIX]{x1D6FC})=\frac{(3-3\unicode[STIX]{x1D6FC}+2d)(1+\unicode[STIX]{x1D6FC})}{4d}\left[1-\frac{1}{64}c^{\ast }(\unicode[STIX]{x1D6FC})\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D708}_{2}^{\ast }(\unicode[STIX]{x1D6FC})=\frac{1+\unicode[STIX]{x1D6FC}}{d-1}\left[\frac{d-1}{2}+\frac{3(d+8)(1-\unicode[STIX]{x1D6FC})}{16}+\frac{4+5d-3(4-d)\unicode[STIX]{x1D6FC}}{1024}c^{\ast }(\unicode[STIX]{x1D6FC})\right], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle c^{\ast }(\unicode[STIX]{x1D6FC})=\frac{32(1-\unicode[STIX]{x1D6FC})(1-2\unicode[STIX]{x1D6FC}^{2})}{9+24d+(8d-41)\unicode[STIX]{x1D6FC}+30\unicode[STIX]{x1D6FC}^{2}(1-\unicode[STIX]{x1D6FC})}, & \displaystyle\end{eqnarray}$$
$d=2$
, since we deal with disks.2.2 The heated granular gas: convective-base-state
First, we revisit the general argument which states that a hydrostatic state is impossible when a temperature gradient is assumed in the horizontal (transverse to gravity) direction, as it is in the case of dissipative lateral walls (DLW) Pontuale et al. (Reference Pontuale, Gnoli, Vega Reyes and Puglisi2016), i.e. when
$T=T(x,y)$
. Momentum balance in (2.2), supplemented by the equation of state for the ideal gas, in the absence of macroscopic flow (hydrostatic) states that
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{x}p=\unicode[STIX]{x2202}_{x}(n(x,y)T(x,y))=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{y}p=\unicode[STIX]{x2202}_{y}(n(x,y)T(x,y))=-mgn(x,y). & \displaystyle\end{eqnarray}$$
$p(x,y)\equiv p(y)$
, which, used in the second equation, sets
$n(x,y)\equiv n(y)$
and, returning this into the first equation above, we obtain
$T(x,y)\equiv T(y)$
, i.e.
$\unicode[STIX]{x2202}T/\unicode[STIX]{x2202}x=0$
. This is in contradiction with the horizontal temperature gradient assumed above. Therefore, the simple hydrostatic system of equations is not compatible with the condition
$\unicode[STIX]{x2202}T/\unicode[STIX]{x2202}x\neq 0$
, originated by the energy sinks at the side walls. We must conclude that a hydrostatic state in the presence of DLW and gravity is not possible. Thus, since
$\boldsymbol{u}\neq 0$
, a flow must always be present, even at infinitesimal values of the Rayleigh number, i.e. there is no hydrostatic solution even if
$Ra\rightarrow 0$
. In the terminology of the previous bibliography on molecular fluids enclosed by dissipative sidewalls, a smooth transition occurs so that the concept of a critical Rayleigh number is no longer tenable (Hall & Walton Reference Hall and Walton1977).
3 Extended Boussinesq-like approximation for a granular gas
However, our previous analysis does not explain why a convection in the form of that observed in the experiments appears. Indeed, our analysis only implies that there is never a hydrostatic solution: it does not necessarily lead to a flow with one convection cell next to each inelastic wall, nor to a convection-free region for points sufficiently far away from the side walls, as seen in the experiments (Wildman et al. Reference Wildman, Huntley and Parker2001a ,Reference Wildman, Huntley and Parker b ; Risso et al. Reference Risso, Soto, Godoy and Cordero2005; Eshuis et al. Reference Eshuis, van der Meer, Alam, van Gerner, van der Weele and Lohse2010; Windows-Yule et al. Reference Windows-Yule, Rivas and Parker2013; Pontuale et al. Reference Pontuale, Gnoli, Vega Reyes and Puglisi2016). Furthermore, we also need to discard it if other properties, other than the sidewalls’ inelasticity, that are present in the experimental system (such as particle–bottom plate friction) are important for the appearance of dissipative lateral walls convection (DLWC) in the granular gas.
Therefore, we need to analyse in § 3 the minimal system of differential equations that derives from the general balance equations and that is able to reproduce a convection with the characteristics of that observed in the experiments. Once it is numerically solved, we will be able to describe in detail the main features of the non-hydrostatic base state in our system.
According to both experiments and computer simulations (molecular dynamics (MD)) results, this convection would only show one cell per inert wall, independent of thermal gradients’ strength and system size (Pontuale et al. Reference Pontuale, Gnoli, Vega Reyes and Puglisi2016). The flow in the bulk of the fluid for wide systems appears to be zero or negligible. This is analogous to the result for molecular fluids where the sidewalls’ effects are important (Hall & Walton Reference Hall and Walton1977). Moreover, the presence of sidewalls introduces two important effects that differ in origin. The first arises from finite-size effects alone and it shows up in a lower critical Rayleigh number for Bénard’s convection, even if the lateral walls are perfectly insulating (i.e. even if they do not convey a lateral energy flux). It is impossible to escape this effect both in molecular (as described in the work by Hall & Walton Reference Hall and Walton1977) and granular fluids when enclosed by lateral walls. The second effect comes properly from energy dissipation at the lateral walls and as we saw produces a non-hydrostatic steady state by default.
We are specifically interested in this second effect that arises only with DLW and leave for future work the study of the first effect (that should eventually lead to consider more appropriate theoretical criteria when comparing with experimental results for the classical bulk-granular convection). Therefore, it is our aim to elucidate the minimal theoretical framework that is able to take into account the important experimental evidence of the DLWC in granular gases.
For our theoretical description, let us use the following reference units: particle mass
$m$
for mass, particle diameter
$\unicode[STIX]{x1D70E}$
for length, thermal velocity at the base
$v_{0}=(4\unicode[STIX]{x03C0}T_{0}/m)^{1/2}$
for velocity, pressure at the base
$p_{0}=n_{0}T_{0}$
, and
$\unicode[STIX]{x1D70E}/v_{0}$
for time, where
$T_{0},n_{0}$
are arbitrary values of temperature and particle density, respectively.
A common situation for thermal convection is that all density derivatives are negligible except for the spatial dependence of density that is coming from gravity, which appears in the momentum balance equation. This happens when the variation of mechanical energy is small compared to the variation of thermal energy (Gray & Giorgini Reference Gray and Giorgini1976) and leads to the Boussinesq equations (Busse Reference Busse1978; Chandrasekhar Reference Chandrasekhar1981). This is always true in our system if the reduced gravitational acceleration fulfils
$g\unicode[STIX]{x1D70E}/v_{0}^{2}\ll 1$
. Thus, we restrict our analysis to small values of
$g$
. Taking this into account in the mass balance equation in (2.2) immediately yields
$\unicode[STIX]{x2202}u_{x}/\unicode[STIX]{x2202}x+\unicode[STIX]{x2202}u_{y}/\unicode[STIX]{x2202}y=0$
. Moreover, for weak convection (as it is the case of our experimental results), we can neglect the advection (nonlinear) terms that emerge in the balance equations (2.2) (Busse Reference Busse1978).
Incorporating these approaches into the other balance equations in (2.2) and for our system geometry (see figure 1) and with our reference units, we get the following dimensionless Boussinesq equations for weak convection in the granular gas
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}^{\ast }(\unicode[STIX]{x1D6FC})\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}\left[\sqrt{T}\left(\frac{\unicode[STIX]{x2202}u_{x}}{\unicode[STIX]{x2202}y}+\frac{\unicode[STIX]{x2202}u_{y}}{\unicode[STIX]{x2202}x}\right)\right]+2\unicode[STIX]{x1D702}^{\ast }(\unicode[STIX]{x1D6FC})\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[\sqrt{T}\frac{\unicode[STIX]{x2202}u_{x}}{\unicode[STIX]{x2202}x}\right]-\frac{\unicode[STIX]{x2202}(nT)}{\unicode[STIX]{x2202}x}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}^{\ast }(\unicode[STIX]{x1D6FC})\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left[\sqrt{T}\left(\frac{\unicode[STIX]{x2202}u_{x}}{\unicode[STIX]{x2202}y}+\frac{\unicode[STIX]{x2202}u_{y}}{\unicode[STIX]{x2202}x}\right)\right]+2\unicode[STIX]{x1D702}^{\ast }(\unicode[STIX]{x1D6FC})\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}\left[\sqrt{T}\frac{\unicode[STIX]{x2202}u_{y}}{\unicode[STIX]{x2202}y}\right]-\frac{\unicode[STIX]{x2202}(nT)}{\unicode[STIX]{x2202}y}-ng^{\ast }=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle n^{2}T^{3/2}\unicode[STIX]{x1D701}^{\ast }(\unicode[STIX]{x1D6FC})=\frac{\unicode[STIX]{x1D705}^{\ast }(\unicode[STIX]{x1D6FC})}{\unicode[STIX]{x03C0}}\left[\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(\sqrt{T}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}x}\right)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}\left(\sqrt{T}\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}y}\right)\right], & \displaystyle\end{eqnarray}$$
with
$g^{\ast }=4\unicode[STIX]{x03C0}g$
. In fact, our Boussinesq-extended approximation includes additional terms with respect to the classical Boussinesq approximation used for thermal convection, since we do not neglect the temperature dependence of the transport coefficients, which results in
$\sqrt{T}$
factors inside the bracket terms in equations (3.1)–(3.3). In a previous work, we noticed that these temperature factors are relevant for important properties of the steady profiles of the granular temperature, such as the curvature (Vega Reyes & Urbach Reference Vega Reyes and Urbach2009). Note that we also keep the density derivatives in (3.1)–(3.2), since the granular gas is highly compressible. For this reason, we do not strictly consider density to be constant in the mass balance equation; we neglect density variations along the flow-field lines instead, while keeping density variations in the momentum balance equation. The corresponding dimensionless forms of the boundary conditions would be:
$T(y=+h/2)=T_{+}/(mv_{0}^{2}/2)$
, (
$[\unicode[STIX]{x2202}T/\unicode[STIX]{x2202}y]_{y=+h/2}=0$
for an open on top system),
$T(y=-h/2)=1/\sqrt{4\unicode[STIX]{x03C0}}$
,
$\boldsymbol{u}(x=\pm L/2)=\boldsymbol{u}(y=\pm h/2)=\mathbf{0}$
,
$p(y=-h/2)=1$
, plus for the lateral walls dissipation condition
$q_{x}(x=\pm L/2)={\mathcal{A}}(\unicode[STIX]{x1D6FC}_{w})[pT^{3/2}]_{x=\pm L/2}$
.
3.1 Numerical solution and comparison with experiments
DLW convection for systems with a top wall and different widths: (a)
$L=h$
, (b)
$L=3h$
and (c)
$L=5h$
. In our dimensionless units: system height
$h=170$
. In each plot, the thickest stream lines correspond to
$u_{0}=(u_{x}^{2}+u_{y}^{2})^{1/2}=3.11\times 10^{-3}$
(in our dimensionless units). The other relevant parameters in this figure are:
$T_{0}=3T_{+}$
and
$g=0.002\,g_{0}$
,
$\unicode[STIX]{x1D6FC}=0.9$
(for both particle–particle and wall–particle collisions).
To numerically solve equations (3.1)–(3.3), we used the finite-volume method. For this, we wrote a code using the SIMPLE algorithm (to avoid numerical decoupling of the pressure field Ferziger & Perić Reference Ferziger and Perić2002) and the FiPy differential equation package with the PySparse solver (Guyer, Wheeler & Warren Reference Guyer, Wheeler and Warren2009). We have seen (see figure 3) that the agreement with experiment and MD simulations is qualitatively very good. In addition, all major properties of the flow are reproduced in the numerical solution obtained from the Boussinesq approximation.
Figure 2, that correspond to systems provided with a top wall, clearly demonstrates that wider systems do not display a greater number of convection cells. The number of cells always remains one per dissipative wall. We also note that the flow is upward in the outer part of the cells (toward the system centre) and downward next to the lateral walls. It is interesting to note also that, according to numerical results, the convection speed
$u_{0}=(u_{x}^{2}+u_{y}^{2})^{1/2}$
reaches roughly the same maximum value when varying the system thickness. Furthermore, if we define the cell size as the horizontal distance between the upward and downward streams (see figure 2) points with maximum convection speed, then we see that the size of the cells remains constant. The exception is for systems thinner than twice the cell size, in which case the cells squeeze each other (figure 2
a). All of these results show the peculiarities of the DLW convection with respect to Rayleigh–Bénard convection and coincide with the experimental behaviour detected previously (Pontuale et al.
Reference Pontuale, Gnoli, Vega Reyes and Puglisi2016).
In figures 3 and 4, that correspond to open systems (i.e. without a top wall), we see a comparison between theory and experimental results, for
$g=0.016\;g_{0}$
(with
$g_{0}=9.8~\text{m}~\text{s}^{-2}$
) in the cases of
$N=300$
and
$1000$
particles in the experimental set-up, respectively. In all cases, we may consider the system as dilute, since the local packing fraction
$\unicode[STIX]{x1D708}=n\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}^{2}/4$
is never greater than approximately
$\unicode[STIX]{x1D708}=0.5\times 10^{-2}$
. As we can see, the agreement is good for the flow field, and more qualitative for the temperature and density fields. In the case of figure 3, there is some disagreement between theory and experiment for both temperature and density, which may be partially explained by the fact that the Knudsen number is not small (
$Kn\sim 10^{-1}$
), and therefore, there might be non-Newtonian effects in the experiments that are not taken into account in our theory. However, it is clear in both cases that the cold fluid regions next to the lateral walls are adjacent to the convective cell centres. We also check that when dissipation at the lateral walls is switched off, no DLW convection appears out of the theoretical solution. In this way, we may conclude that the convection mechanism appears as a consequence of the combined action of two perpendicular gradients: the density (and thermal) gradient due to the action of gravity and the horizontal gradient due to energy dissipation by the lateral walls.
Let us point out that our theory does not take into account all of the ingredients and details that are present in previous experiments on DLWC in granular systems (Wildman et al. Reference Wildman, Huntley and Parker2001a ; Windows-Yule et al. Reference Windows-Yule, Rivas and Parker2013; Pontuale et al. Reference Pontuale, Gnoli, Vega Reyes and Puglisi2016), such as plate–particle friction and/or sliding or dynamical effects derived from particle sphericity (just to give two examples), etc. Moreover, as noted in previous works, there is a tendency for volume convection disappearance for not-so-small Knudsen numbers (Ansari & Alam Reference Ansari and Alam2016; Pontuale et al. Reference Pontuale, Gnoli, Vega Reyes and Puglisi2016). For all of this, a comparison with our previous experimental results (Pontuale et al. Reference Pontuale, Gnoli, Vega Reyes and Puglisi2016) and the experimental results by others (Wildman et al. Reference Wildman, Huntley and Parker2001a ; Eshuis et al. Reference Eshuis, van der Meer, Alam, van Gerner, van der Weele and Lohse2010; Windows-Yule et al. Reference Windows-Yule, Rivas and Parker2013) should be regarded as qualitative, not quantitative. The advantage of this procedure is that, because it is reduced to the essentials, we are finally able to identify the key ingredients that produce the DLWC in the granular gas.
Fields from experiments and theory (left and right panels, respectively, for each pair of panels for the corresponding field), for a system without a top wall (open system). Length unit: particles diameter (
$\unicode[STIX]{x1D70E}=1~\text{mm}$
). Here
$g=0.016\times 9.8~\text{m}~\text{s}^{-2}$
. (a,b) Flow field (
$\boldsymbol{u}=\mathbf{0}$
). (c,d) The corresponding temperature
$T/m$
for the experiment (c) and theory (d). (e,f) Packing fraction fields
$\unicode[STIX]{x1D708}=n\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}^{2}/4$
for the experiment (e) and theory (f). Black stands for lower and red for higher field values. Here
$0<T\leqslant 0.2mv_{0}^{2}$
,
$(v_{0}=370~\text{mm}~\text{s}^{-1})$
;
$0.02\times 10^{-2}<\unicode[STIX]{x1D708}\leqslant 0.15\times 10^{-2}$
, mean packing fraction
$\unicode[STIX]{x1D708}=0.049\times 10^{-2}$
. Density colour bars are in percentage units. Experiments: with
$f=45~\text{Hz}$
,
$A=1.85~\text{mm}$
,
$N=300$
(total number of particles). Theory: coefficient of restitution
$\unicode[STIX]{x1D6FC}=0.9$
(both for particle–particle and wall–particle collisions).
Fields from experiments and theory (left and right panels, for each pair of panels for the corresponding field, respectively) for a system without a top wall (open system). Length unit: particles diameter (
$\unicode[STIX]{x1D70E}=1~\text{mm}$
). Here
$g=0.016\times 9.8~\text{m}~\text{s}^{-2}$
. (a,b) Flow field (
$\boldsymbol{u}=\mathbf{0}$
). (c,d) The corresponding temperature
$T/m$
for experiment and theory (d). (e,f) Packing fraction
$\unicode[STIX]{x1D708}=n\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}^{2}/4$
fields for experiment (e) and theory (f). Black denotes lower and red denotes higher field values. Here
$0<T\leqslant 0.2mv_{0}^{2}$
,
$(v_{0}=370~\text{mm}~\text{s}^{-1})$
;
$0.02\times 10^{-2}<\unicode[STIX]{x1D708}\leqslant 0.46\times 10^{-2}$
, mean packing fraction
$\unicode[STIX]{x1D708}=0.15\times 10^{-2}$
. Density colour bars are in percentage units. Experiments: with
$f=45~\text{Hz}$
,
$A=1.85~\text{mm}$
,
$N=1000$
(total number of particles). Theory: coefficient of restitution
$\unicode[STIX]{x1D6FC}=0.9$
(both for particle–particle and wall–particle collisions).
Furthermore, theoretical procedures (simulations, hydrodynamics) allow us to separate the two sources of dissipation and make ‘ideal’ assumptions. For instance, we can switch off internal dissipation (and reproduce the molecular fluid limit case) while retaining dissipation at the walls.
4 Conclusions
In this work we have discussed the theory framework for a previously observed experimental phenomenon of granular convection (Wildman et al.
Reference Wildman, Huntley and Parker2001a
,Reference Wildman, Huntley and Parker
b
; Risso et al.
Reference Risso, Soto, Godoy and Cordero2005; Eshuis et al.
Reference Eshuis, van der Meer, Alam, van Gerner, van der Weele and Lohse2010; Windows-Yule et al.
Reference Windows-Yule, Rivas and Parker2013; Pontuale et al.
Reference Pontuale, Gnoli, Vega Reyes and Puglisi2016). This convection appears automatically, i.e. occurs at arbitrary
$Ra$
and, in close analogy to the convection in molecular fluids induced by cold sidewalls (Hall & Walton Reference Hall and Walton1977), we conclude that the new granular convection is also induced by dissipative vertical sidewalls and a gravity field. We have denoted it as a DLWC. We have built a granular hydrodynamic theoretical framework that explains the physical origin of this convection, in the context of the Boussinesq equations for the granular gas. To the best of the authors’ knowledge, this is the first time that a Boussinesq-type approach has been used for granular convection. We have found that our theory also explains the main features of the experimental observations. That is, the DLW convection displays, in all cases, only one convection cell per dissipative wall, the width of this wall being increased when convection intensity increases. Moreover, the DLW convection intensity is enhanced by increasing gravity acceleration, wall dissipation, or both. Conversely, it is decreased by increasing bottom wall temperature (at fixed gravity and wall dissipation).
Note that sidewalls are inherently inelastic in granular fluid experimental systems. Thus, the DLWC is always present at the experimental level and our results imply that the classical volume thermal convection in granular gases does not appear in experimental systems out of a hydrostatic state. Instead, it develops as a secondary instability out of the DLWC state. Therefore, more theory work would be needed in general to correctly describe the experimental instability criteria for the volume thermal convection in granular gases. This also implies that the accuracy of previous hydrodynamic theory for the volume thermal convection can be seemingly improved if the prior existence of the DLW is taken into account.
Finally, our present work constitutes further strong evidence that steady granular flows can be described correctly by a standard hydrodynamic theory (Puglisi Reference Puglisi2015) (see also the recent work by Vega Reyes & Lasanta (Reference Vega Reyes and Lasanta2017), with results in the same line).
Acknowledgements
F.V.R. acknowledges support from the Salvador Madariaga program and grant numbers PRX14/00137, FIS2016-76359-P (both from the Spanish Government) and GR15104 (Regional Government of Extremadura, Spain). Use of computing facilities from the Extremadura Research Centre for Advanced Technologies (CETA-CIEMAT) partially financed by the ERDF is also acknowledged.
