Vadasz Number Effects on Convection in a Horizontal Porous Layer Subjected to Internal Heat Generation and G-Jitter

Abstract

Flow and heat transfer in a horizontal porous layer subjected to internal heat generation and g-jitter is considered for the Dirichlet thermal boundary condition. A linear stability analysis is used to determine the convection threshold in terms of the critical Rayleigh number. For the low amplitude, high frequency approximation, the results show that vibration has a stabilizing effect on the onset of convection when the porous layer is heated from below. When the porous layer is cooled from below and heated from above, the vibration has a destabilizing effect in the presence of internal heat generation. It is also demonstrated that when the top and bottoms walls are cooled and rigid/impermeable, the critical Rayleigh number is infinitely large and conduction is the only possible mode of heat transfer. The impact of increasing the Vadasz number is to stabilize the convection, in addition to reducing the transition point from synchronous to subharmonic solutions.

1. Introduction

Gravity driven natural convection in fluid saturated porous media has been widely studied for various configurations for proposed engineering applications. Applications include inter alia rotating machinery, chemical and materials processing, food industry and possibly nuclear vessel technology. Pioneering past works, where the basic temperature gradient is constant, have been investigated by several authors. Early works by [1,2,3,4] for non-rotating porous media provides an excellent reference base for readers for constant basic temperature gradient. These early works set the foundational basis for researchers in heat and mass transfer in porous media. Later important works involving rotating effects in porous media were presented in [5,6,7,8,9,10,11], with [9] providing important results involving Coriolis effects for the case when the centrifugal effects are negligible. The coefficient of the time derivative in the momentum equation in [9] was proposed by [12] to be called the Vadasz number. Typically the Vadasz number is very large in porous media applications and for this reason in most formulations, the impact of the time derivative in the momentum equation is negligible and therefore omitted. The work of [9] considers scenarios where the Vadasz number could be small, thus motivating for the retention of the time derivative in the momentum equation. The results presented in [9] show that the oscillatory mode of convection is possible when the Vadasz number is small to moderate. The work presented in [13,14,15,16] involves the effects of g-jitter on the stability of convection in porous media. In the studies [13,14], the author clearly demonstrates that the transition from synchronous to subharmonic solutions is characterized by a spike in the curve of the synchronous mode just prior to the transition.

The focus of the current work is to consider the effects of vibration (g-jitter) for a porous layer with internal heat generation. In principle this case corresponds to one where the basic temperature is non-constant but rather a function of the vertical z- co-ordinate. Researchers are referred to earlier works for flow and heat transfer with internal heat generation with constant gravity [17,18] and also with variable gravity [19]. The work presented in [17] adopts the Dirichlet thermal boundary condition and the results are quite relevant to the current study. Therefore outputs of the current work will be compared to [17].

Early work on vibration’s effects on porous media with internal heat generation are presented in [20,21]. Later analyses of vibration effects in a porous layer with through flow and internal heat generation are presented in [22], which is an extension of previous work by this author [23,24]. Other works including internal heat generation in porous media with anisotropy are presented in [25,26]. Research specifically involving the effects of gravity modulation and internal heat generation in porous media are presented in [25,26,27,28]. Readers are also referred to two important books [29,30] for further review/reading.

The current study has applications in nuclear reactor technology where possible seismic effects may be modelled in terms of the g-jitter term. As an example, pebble bed modular reactors used in gas cooled reactor plants with helium turbines could serve as a specific example of an engineering application directly aligned with the current model being analysed. The gas-cooled reactor is typically used in closed systems that are made up of a recuperator, helium turbine and inventory control/storage systems. When this type of plant is operated, the circulating helium cools the porous pebble bed reactor and in doing so heats up to a maximum of about 1100 °C. The heated helium gas then passes to a helium turbine to produce power. Studies around porous media subjected to internal heat generation and vibration effects for high temperature reactors should therefore largely be driven by the safety aspects surrounding this technology [17]. Given that the future energy generation mix will focus largely on low carbon emissions technologies, with significantly reduced water usage, there is a renewed interest in pebble bed and molten salt reactor technology, hence the impetus for the current work. Although the analysis in this paper uses a gas cooled porous reactor as an example, other examples for application of this work include systems involving heat extraction or storage (in porous media) using molten salts. These storage systems are largely found in the concentrating solar plant (CSP), and more specifically central receiver systems using molten salts as the working fluid. The CSP plant is typically made up of the heliostat field, the receiver tower where the solar rays are focused to heat up the salt and the steam generator that exchanges heat from the hot salt loop to the bottoming Rankine steam cycle. Some of the heated salt is also directed to the storage tanks containing porous media. In CSPs, maximizing hours of storage is critical to increase the plant capacity factor. The size of storage (contributing to capacity factor) and ultimately the solar multiple for the plant are traded off by minimizing the levelised cost for electricity. Therefore a better understanding of storage systems and modelling thereof, with a view to optimizing them to operate efficiently becomes critical, due to the impact on levelised cost of electricity. The storage tank in a typical solar plant is made up of a large circular packed porous bed that has permeable upper and lower boundaries. The porous bed in the tank is saturated by a high freezing point ammonium salt. Further research and development in this area will make the plant more efficient from a heat storage perspective. Also, seismic effects and an understanding of these effects on the transfer of heat under these conditions become important as a design/containment consideration to protect the environment.

The current study will determine the basic conduction flow and temperature profiles and then solve the characteristic Rayleigh number definition using the linear stability theory. Following [13,31], the equations that result will be cast into the canonical form of the Mathieu equations and thereafter the critical Rayleigh number in terms of the Vadasz number will be determined for the case of low-amplitude high-frequency vibration. The data will then be compared to [17], and the case of a porous layer subjected to vibration without internal heat generation. This study is distinctly different from [20,21,22] and will not resort to the use of a numerical solver, as the intent is to recover the design equations as a function of the key contributing variables. Based on the review performed, the author is unable to locate any paper that has adopted the current methodology or for that matter has produced the results presented.

2. Problem Formulation

Figure 1 shows a horizontal fluid saturated porous layer subjected to vibration that forms the basis for the current work. The porous layer is sandwiched between two rigid and impermeable horizontal plates, (a distance ${\mathrm{H}}^{*}$ apart), subjected to gravity ${\mathrm{g}}^{*}$, and vibration ${\mathrm{b}}^{*}{\mathrm{\omega }}^{*2}\sin \left({\mathrm{\omega }}^{*}{\mathrm{t}}^{*}\right)$. In Figure 1, the x* refers to the horizontal distances and z* refers the vertical distances. At ${\mathrm{z}}^{*}=0$: $\partial {\mathrm{T}}^{*}/\partial {\mathrm{z}}^{*}=0$ and at ${\mathrm{z}}^{*}=1$:${\mathrm{T}}^{*}={\mathrm{T}}_{\mathrm{C}}^{*}$ which represents the adiabatic lower wall and perfectly conducting upper wall. The Boussinesq approximation is applied to model the effects of the density variations and results in the following system of dimensional equations for continuity, momentum and energy:

$${\nabla }^{*}\cdot {\mathbf{V}}^{*}=0,$$

(1)

$$\left(\frac{{\mathrm{\rho }}_{\mathrm{o}}^{*}}{{\mathrm{\varphi }}^{*}}\frac{\partial }{\partial {\mathrm{t}}^{*}}+1\right){\mathbf{V}}^{*}=\frac{{\mathrm{\chi }}^{*}}{{\mathrm{\mu }}^{*}}\left[-{\nabla }^{*}{\mathrm{p}}^{*}+{\mathrm{\rho }}^{*}\left({\mathrm{g}}^{*}+{\mathrm{b}}^{*}{\mathrm{\omega }}^{*2}\sin \left({\mathrm{\omega }}^{*}{\mathrm{t}}^{*}\right)\right)\right]{\widehat{\mathbf{e}}}_{\mathbf{z}},$$

(2)

$$\frac{\partial {\mathrm{T}}^{*}}{\partial {\mathrm{t}}^{*}}+{\mathbf{V}}^{*}\cdot {\nabla }^{*}{\mathrm{T}}^{*}={\mathrm{\kappa }}^{*}\left({\nabla }^{*2}{\mathrm{T}}^{*}+\frac{{\mathrm{q}}^{*}}{{\mathrm{k}}_{\mathrm{m}}^{*}}\right).$$

(3)

In the system (1–3), the fluid velocity is ${\mathbf{V}}^{*}$, the pressure is ${\mathrm{p}}^{*}$ and the density is ${\mathrm{\rho }}^{*}$. In Equations (1)–(3), the temperature is ${\mathrm{T}}^{*}$, the permeability of the porous medium is ${\mathrm{\chi }}^{*}$, the dynamic viscosity is ${\mathrm{\mu }}^{*}$ and the fluid thermal conductivity is ${\mathrm{k}}_{\mathrm{m}}^{*}$. The governing equations may be non-dimensionalised using the scaling variables ${\mathrm{\kappa }}^{*}/{\mathrm{L}}^{*}$, ${\mathrm{\mu }}^{*}{\mathrm{\kappa }}^{*}/{\mathrm{k}}_0^{*}$ and $\mathsf{\Delta }{\mathrm{T}}^{*}={\mathrm{T}}_{\mathrm{H}}^{*}-{\mathrm{T}}_{\mathrm{C}}^{*}$ for the filtration velocity components $\left({\mathrm{u}}^{*},{\mathrm{v}}^{*},{\mathrm{w}}^{*}\right)$, reduced pressure $\left({\mathrm{p}}^{*}\right)$, and the temperature variations ${\mathrm{T}}^{*}-{\mathrm{T}}_{\mathrm{C}}^{*}$. In these scalings, ${\mathrm{\kappa }}^{*}$ is the thermal diffusivity and ${\mathrm{k}}_0^{*}$ is a characteristic permeability associated with the porous medium. The height of the porous medium ${\mathrm{H}}^{*}$ is used to scale all spatial lengths according to the relations, $\mathrm{x}={\mathrm{x}}^{*}/{\mathrm{H}}^{*}$, $\mathrm{y}={\mathrm{y}}^{*}/{\mathrm{H}}^{*}$ and $\mathrm{z}={\mathrm{z}}^{*}/{\mathrm{H}}^{*}$. Applying the scaling factors to Equations (1)–(3) yields the following for constant permeability, in other words, $\mathrm{\chi }=1$:

$$\nabla \cdot \mathbf{V}=0,$$

(4)

$$\left(\frac{1}{\mathrm{V}\mathrm{a}}\frac{\partial }{\partial \mathrm{t}}+1\right)\mathbf{V}=\left[-\nabla \mathrm{p}+\mathrm{R}{\mathrm{a}}_{\mathrm{g}}\left(1+\mathrm{\delta }\sin \left(\mathrm{\omega }\mathrm{t}\right)\right)\right]\mathrm{T}{\widehat{\mathrm{e}}}_{\mathrm{z}},$$

(5)

$$\frac{\partial \mathrm{T}}{\partial \mathrm{t}}+\mathbf{V}\cdot \nabla \mathrm{T}={\nabla }^2\mathrm{T}+1,$$

(6)

The key non-dimensional parameters that emanate from the rescaling of Equations (1)–(3) are the gravitational Rayleigh number, $\mathrm{R}{\mathrm{a}}_{\mathrm{g}}=\left({\mathrm{\rho }}_0^{*}{\mathrm{\beta }}^{*}{\mathrm{g}}^{*}{\mathrm{H}}^{*}{\mathrm{k}}_0^{*}/\left({\mathrm{\mu }}^{*}{\mathrm{\kappa }}^{*}\right)\right)\left({\mathrm{q}}^{*}{\mathrm{H}}^{*2}/{\mathrm{k}}_{\mathrm{m}}^{*}\right)$, the vibration amplitude $\mathrm{\delta }=\mathrm{\eta }\mathrm{F}\mathrm{r}{\mathrm{\omega }}^2$ where $\mathrm{\eta }={\mathrm{b}}^{*}/{\mathrm{H}}^{*}$ and the Froude number $\mathrm{F}\mathrm{r}={\mathrm{\kappa }}^{*2}/{\mathrm{g}}^{*}{\mathrm{H}}^{*3}$. The Rayleigh number definition in this paper is somewhat different from that defined in [17] and this has to do largely with the variables chosen when non-dimensionalizing the governing equations. In particular readers will observe from [17], that the scaling of the temperature is different from the current work. It will be demonstrated later how this impacts the numerical values of the Rayleigh number when compared to [17]. In Equation (5), Va is the Vadasz number defined as $\mathrm{V}\mathrm{a}={\mathrm{\varphi }}^{*}\mathrm{Pr}/\mathrm{D}\mathrm{a}$, (where Pr is the Prandtl number and Da is the Darcy number) and the symbols V, T and ${\mathrm{p}}_{\mathrm{r}}$ represent the dimensionless filtration velocity vector, temperature and reduced pressure, respectively, whilst ${\widehat{\mathbf{e}}}_{\mathbf{z}}$, is the unit vector in the z- direction. In the previous work [9], the author has put forward a motivation for specific cases when the Vadasz number is small and can be retained. In the instance of liquid metals this could be the case. In the current study when considering a high temperature reactor for the proof of concept geometry, an average porosity of ${\mathrm{\varphi }}^{*}\approx 0.4$, a Prandtl number $\mathrm{Pr}=0.7$, characteristic permeability ${\mathrm{k}}_0^{*}\approx 0.704$ and calculated Darcy number $\mathrm{D}\mathrm{a}\approx 0.155$ yields a Vadasz number $\mathrm{V}\mathrm{a}=1.81$. This implies that for systems involving gas reactors and porous media, there are instances in which the Vadasz number is close to unity and can be retained in the momentum equation. The solutions for the basic temperature and flow field is given as ${\mathrm{T}}_{\mathrm{B}}=\mathrm{B}+\mathrm{A}\mathrm{z}-1/2{\mathrm{z}}^2$ and ${\mathrm{V}}_{\mathrm{B}}=0$, where A and B are constants that may be determined, based on the imposed boundary for the various cases shown in

Table 1. Configurations for upper and lower wall heating with internal heat generation.

Case	Boundary Condition	A	Description
1	$z=0:\frac{d{T}_B}{dz}=0$$z=1:T_B=0$	0	Adiabatic bottom wall/Perfectly conducting top wall
2	$z=0:T_B=0$$z=1:T_B=0$	$\frac{1}{2}$	Perfectly conducting top and bottom walls
3	$z=0:T_B=1$$z=1:T_B=0$	$-\frac{1}{2}$	Perfectly conducting top and bottom walls
4	$z=0:T_B=0$$z=1:T_B=1$	$\frac{3}{2}$	Perfectly conducting top and bottom walls
5	$z=0:T_B=0$$z=1:\frac{d{T}_B}{dz}=0$	1	Adiabatic top wall/Perfectly conducting bottom wall

. Since we will consider the derivative of the basic temperature in the energy equation, the value for the constant B is not required, and therefore will not be presented in Table 1.

In the current study we will only consider the Dirichlet thermal boundary conditions, in other words, cases 2–4. The Neumann thermal boundary conditions, in other words, cases 1 and 5 will be dealt with separately in another study. The results for the thermal boundary conditions will be compared to the work of [17] and also referenced to the case of Benard convection with g-jitter and no internal heat generation.

In providing a solution to the system of equations, it is convenient to apply the curl operator ($\nabla \times$) twice on Equation (5), and using Equation (4), consider only the vertical z-component of the result to yield,

$$\left(\frac{1}{\mathrm{V}\mathrm{a}}\frac{\partial }{\partial \mathrm{t}}+1\right){\nabla }^2\mathrm{w}-\mathrm{R}{\mathrm{a}}_{\mathrm{g}}\left(1+\mathrm{\delta }\sin \left(\mathrm{\omega }\mathrm{t}\right)\right){\nabla }_{\mathrm{H}}^2\mathrm{T}=0,$$

(7)

where ${\nabla }_{\mathrm{H}}^2={\partial }^2/\partial {\mathrm{x}}^2+{\partial }^2/\partial {\mathrm{y}}^2$.

3. Linear Stability Analysis

Assuming small perturbations around the basic solution of the form $\mathrm{T}={\mathrm{T}}_{\mathrm{B}}+{\mathrm{T}}^{\prime }$ and $\mathrm{w}={\mathrm{w}}_{\mathrm{B}}+{\mathrm{w}}^{\prime }$ and linearizing Equations (4)–(7) yield the following equations,

$$\left(\frac{1}{\mathrm{V}\mathrm{a}}\frac{\partial }{\partial \mathrm{t}}+1\right){\nabla }^2{\mathrm{w}}^{\prime }-\mathrm{R}{\mathrm{a}}_{\mathrm{g}}\left(1+\mathrm{\delta }\sin \left(\mathrm{\omega }\mathrm{t}\right)\right){\nabla }_{\mathrm{H}}^2{\mathrm{T}}^{\prime }=0,$$

(8)

$$\frac{\mathrm{d}{\mathrm{T}}_{\mathrm{B}}}{\mathrm{d}\mathrm{z}}{\mathrm{w}}^{\prime }=\left(\frac{\partial {\mathrm{T}}^{\prime }}{\partial \mathrm{t}}-{\nabla }^2{\mathrm{T}}^{\prime }\right).$$

(9)

Substituting Equation (9) into Equation (8) and simplifying, yields the following equation for ${\mathrm{T}}^{\prime }$,

$$2\left(\frac{1}{\mathrm{V}\mathrm{a}}\frac{\partial }{\partial \mathrm{t}}+1\right)\left({\nabla }^2-\frac{\partial }{\partial \mathrm{t}}\right){\mathrm{T}}^{\prime }+{\left(\mathrm{A}-\mathrm{z}\right)}^2\left(\frac{1}{\mathrm{V}\mathrm{a}}\frac{\partial }{\partial \mathrm{t}}+1\right){\nabla }^2\left({\nabla }^2-\frac{\partial }{\partial \mathrm{t}}\right){\mathrm{T}}^{\prime }=\mathrm{R}{\mathrm{a}}_{\mathrm{g}}{\left(\mathrm{A}-\mathrm{z}\right)}^3\left(1+\mathrm{\delta }\sin \left(\mathrm{\omega }\mathrm{t}\right)\right){\nabla }_{\mathrm{H}}^2{\mathrm{T}}^{\prime }.$$

(10)

We assume an expansion into normal modes in the x- and y- directions of the form,

$${\mathrm{T}}^{\prime }={\mathrm{e}}^{\mathrm{i}\left({\mathrm{s}}_{\mathrm{x}}\mathrm{x}+{\mathrm{s}}_{\mathrm{y}}\mathrm{y}\right)}{\displaystyle \sum _{\mathrm{k}=1}^{\mathrm{N}}{\mathrm{a}}_{\mathrm{k}}}\left(\mathrm{t}\right)\sin (\mathrm{k}\mathrm{\pi }\mathrm{z}),$$

(11)

which satisfies the boundary conditions, ${\mathrm{T}}^{\prime }=0$ and ${\mathrm{w}}^{\prime }=0$ at $\mathrm{z}=0$ and $\mathrm{z}=1$. In Equation (10) ${\mathrm{s}}^2={\mathrm{s}}_{\mathrm{y}}^2+{\mathrm{s}}_{\mathrm{z}}^2$ and k = 1, 2, 3,…N. Substituting Equation (11) into Equation (10) and then multiplying the result by orthogonal functions and integrating over $\mathrm{z}\in \left[0,1\right]$ yields the following system of equations,

$$\begin{gathered} {\displaystyle \sum _{\mathrm{k}=1}^{\mathrm{N}}\{\left(\frac{1}{\mathrm{V}\mathrm{a}}\frac{\partial }{\partial \mathrm{t}}+1\right)}\left({\mathrm{k}}^2{\mathrm{\pi }}^2+{\mathrm{s}}^2+\frac{\partial }{\partial \mathrm{t}}\right)\left(-1+\left({\mathrm{k}}^2{\mathrm{\pi }}^2+{\mathrm{s}}^2\right)\left(\frac{\mathrm{A}}{2}\left(\mathrm{A}-1\right)+\frac{1}{6}-\frac{1}{{\mathrm{\pi }}^2{(\mathrm{k}+\mathrm{l})}^2}\right)\right)+ \\ +\mathrm{R}{\mathrm{a}}_{\mathrm{g}}{\mathrm{s}}^2\left(1+\mathrm{\delta }\sin \left(\mathrm{\omega }\mathrm{t}\right)\right)\left(\frac{\left(\mathrm{A}-1\right)}{2}\left({\mathrm{A}}^2-\frac{\mathrm{A}}{2}+\frac{1}{2}\right)+\frac{1}{8}+\frac{3}{2}\frac{\left(1-2\mathrm{A}\right)}{{\mathrm{\pi }}^2{\left(\mathrm{l}+\mathrm{k}\right)}^2}\right)\}{\mathrm{\xi }}_{\mathrm{l}\mathrm{k}}{\mathrm{a}}_{\mathrm{k}}+ \\ {\displaystyle \sum _{\mathrm{k}=1}^{\mathrm{N}}\{\left({\mathrm{k}}^2{\mathrm{\pi }}^2+{\mathrm{s}}^2\right)\left(\frac{1}{\mathrm{V}\mathrm{a}}\frac{\partial }{\partial \mathrm{t}}+1\right)}\left({\mathrm{k}}^2{\mathrm{\pi }}^2+{\mathrm{s}}^2+\frac{\partial }{\partial \mathrm{t}}\right)\frac{2\mathrm{k}\mathrm{l}\left(2\mathrm{A}-1\right)}{{\mathrm{\pi }}^2{\left({\mathrm{k}}^2-{\mathrm{l}}^2\right)}^2}+, \\ +\mathrm{R}{\mathrm{a}}_{\mathrm{g}}{\mathrm{s}}^2\left(1+\mathrm{\delta }\sin \left(\mathrm{\omega }\mathrm{t}\right)\right)\left(\frac{6\mathrm{k}\mathrm{l}{\left(\mathrm{A}-1\right)}^2}{{\mathrm{\pi }}^2}-\frac{3\mathrm{k}\mathrm{l}\left(1-2\mathrm{A}\right)}{{\mathrm{\pi }}^2}-\frac{6}{{\mathrm{\pi }}^4}\left(\frac{1}{{\left(\mathrm{l}-\mathrm{k}\right)}^4}-\frac{1}{{\left(\mathrm{l}+\mathrm{k}\right)}^4}\right)\right)\}{\mathrm{\xi }}_{\mathrm{l}+\mathrm{k},2\mathrm{p}-1}{\mathrm{a}}_{\mathrm{k}}=0, \end{gathered}$$

(12)

where ${\mathrm{\xi }}_{\mathrm{l}+\mathrm{k},2\mathrm{p}-1}$ is 1 when $\left(\mathrm{l}+\mathrm{k}\right)$ is odd and zero otherwise. As can be observed, the system shown in Equation (12) is complicated when solved to higher ranks of N. One may develop a numerical solver to solve Equation (12) to the higher ranks, but that is outside the scope of the current work. Even though there may be some inaccuracies associated with considering Equation (12) to rank N = 1, useful information may be drawn from rank N = 1. As one proceeds to higher ranks, one would expect to see a convergence in the critical Rayleigh number. For the current study, considering rank N = 1, Equation (12) may be presented as,

$$\frac{{\mathrm{d}}^2{\mathrm{a}}_1}{\mathrm{d}{\mathrm{t}}^2}+\left(1+\mathrm{V}\mathrm{a}\left({\mathrm{s}}^2+{\mathrm{\pi }}^2\right)\right)\frac{\mathrm{d}{\mathrm{a}}_1}{\mathrm{d}\mathrm{t}}+\mathrm{V}\mathrm{a}{\mathrm{F}}_0{\mathrm{s}}^2\left(\frac{\left({\mathrm{s}}^2+{\mathrm{\pi }}^2\right)}{{\mathrm{F}}_0{\mathrm{s}}^2}+\mathrm{R}{\mathrm{a}}_{\mathrm{g}}\left(1+\mathrm{\delta }\sin \left(\mathrm{\omega }\mathrm{t}\right)\right)\right){\mathrm{a}}_1=0,$$

(13a)

where

$${\mathrm{F}}_0=\frac{4\left(\mathrm{A}-1\right)\left({\mathrm{A}}^2-\frac{1}{2}\mathrm{A}+\frac{1}{2}\right)+1+\frac{3\left(1-2\mathrm{A}\right)}{{\mathrm{\pi }}^2}}{\left(\left(\frac{4}{3}-\frac{2}{{\mathrm{\pi }}^2}\right)+4\mathrm{A}\left(\mathrm{A}-1\right)\right)\left({\mathrm{s}}^2+{\mathrm{\pi }}^2\right)-8},$$

(13b)

$$2\mathrm{p}=1+\mathrm{V}\mathrm{a}\left({\mathrm{s}}^2+{\mathrm{\pi }}^2\right).$$

(13c)

Equation (13a) may be transformed into the Mathieu equation by taking ${\mathrm{a}}_1={\mathrm{e}}^{-\mathrm{\lambda }\mathrm{\tau }}{\mathrm{X}}_1\left(\mathrm{\tau }\right)$ so that the resulting equation may be presented as,

$${{\mathrm{X}}^{″}}_1+\left[\mathrm{M}+2\mathrm{Q}\cos \left(2\mathrm{\tau }\right)\right]{\mathrm{X}}_1=0.$$

(14)

In Equation (14), the coefficients M and 2Q are defined for stationary convection as,

$$\mathrm{M}=\frac{4\mathrm{V}\mathrm{a}{\mathrm{F}}_0{\mathrm{s}}^2}{{\mathrm{\omega }}^2}\left(\mathrm{R}{\mathrm{a}}_{\mathrm{g}}+\frac{{\mathrm{s}}^2+{\mathrm{\pi }}^2}{{\mathrm{F}}_0{\mathrm{s}}^2}\right),$$

(15)

$$2\mathrm{Q}=\frac{4\mathrm{V}\mathrm{a}{\mathrm{F}}_0{\mathrm{s}}^2\mathrm{R}{\mathrm{a}}_{\mathrm{g}}\mathrm{\delta }}{{\mathrm{\omega }}^2}.$$

(16)

In previous studies [14,20], authors have related M and Q in Equations (15)–(16) via an indirect numerical method. Whilst that method recovers both the synchronous and subharmonic modes, it is quite cumbersome. In this paper I propose the following asymptotic expansions for low amplitude and high frequency vibrations in an attempt to link M and Q,

$${\mathrm{X}}_1={\mathrm{X}}_0+\mathrm{Q}{\mathrm{X}}_1+{\mathrm{Q}}^2{\mathrm{X}}_2+\dots ,$$

(17)

$$\mathrm{M}={\mathrm{M}}_0+\mathrm{Q}{\mathrm{M}}_1+{\mathrm{Q}}^2{\mathrm{M}}_2+\dots ,$$

(18)

where Q as defined in Equation (16) assumes small values for high vibration frequencies. The methodology used in the derivation of Equations (19) and (20) involves the use of simple asymptotic expansions, in other words, Equations (17) and (18) that are then applied on the Mathieu equation, Equation (14). The resulting equations are solved to the various orders in Q noting that one needs to apply the solvability condition to the resulting differential equations. The result shown in Equation (20) therefore is the solvability condition. The relation between M and Q are then used in conjunction with their definitions in Equations (15) and (16) to yield the characteristic equation for the Rayleigh number. In the derivation of Equations (19) and (20), the coefficient of exponential growth terms are set to zero and the constant, ${\mathrm{b}}_0>0$, is an integration constant and a real number. The resulting solutions to rank N = 1 and the solvability condition may be represented as,

$${\mathrm{X}}_1={\mathrm{c}}_0\left(1-\frac{\mathrm{Q}}{2}\cos \left(2\mathrm{\tau }\right)\right),$$

(19)

$$\mathrm{M}=-\frac{1}{2}{\mathrm{Q}}^2.$$

(20)

Using Equations (15) and (16) and substituting into Equation (20) yields the following characteristic equation for the Rayleigh number,

$$\frac{1}{2}{\mathrm{F}}_1\mathrm{V}\mathrm{a}{\mathrm{s}}^2{\left(\mathrm{\kappa }\mathrm{F}\mathrm{r}\mathrm{\omega }\right)}^2\mathrm{R}{\mathrm{a}}_{\mathrm{g}}^2-\mathrm{R}{\mathrm{a}}_{\mathrm{g}}+\mathrm{R}{\mathrm{a}}_{\mathrm{g}0}=0$$

(21a)

where

$${\mathrm{F}}_1=\frac{\frac{1}{2}\left(1-\mathrm{A}\right)\left({\mathrm{A}}^2-\frac{1}{2}\mathrm{A}+\frac{1}{2}\right)-\frac{1}{8}-\frac{3}{8{\mathrm{\pi }}^2}\left(1-2\mathrm{A}\right)}{\left(\left(\frac{1}{6}-\frac{1}{4{\mathrm{\pi }}^2}\right)+\frac{1}{2}\mathrm{A}\left(\mathrm{A}-1\right)\right)\left({\mathrm{s}}^2+{\mathrm{\pi }}^2\right)-1},$$

(21b)

$$\mathrm{R}{\mathrm{a}}_{\mathrm{g}0}=\frac{\left({\mathrm{s}}^2+{\mathrm{\pi }}^2\right)\left(\left(\left(\frac{1}{6}-\frac{1}{4{\mathrm{\pi }}^2}\right)+\frac{1}{2}\mathrm{A}\left(\mathrm{A}-1\right)\right)\left({\mathrm{s}}^2+{\mathrm{\pi }}^2\right)-1\right)}{{\mathrm{s}}^2\left(\frac{1}{2}\left(\left(1-\mathrm{A}\right)\left({\mathrm{A}}^2-\frac{1}{2}\mathrm{A}+\frac{1}{2}\right)\right)-\frac{1}{8}-\frac{3}{8{\mathrm{\pi }}^2}\left(1-2\mathrm{A}\right)\right)}.$$

(21c)

From Equation (21a) one may observe that for the cases when the vibration frequency approaches zero, in other words, $\mathrm{\omega }\to 0$, the Rayleigh number approaches the definition, $\mathrm{R}{\mathrm{a}}_{\mathrm{g}}\to \mathrm{R}{\mathrm{a}}_{\mathrm{g}0}$. Equation (21a) is solved for both zero and non-zero vibration frequencies for case 2 to case 4.

4. Results and Discussion

The Vadasz number and the Rayleigh numbers are scaled as follows: $\mathrm{\gamma }=\mathrm{V}\mathrm{a}/{\mathrm{\pi }}^2$, $\widehat{\mathrm{R}}=\mathrm{R}{\mathrm{a}}_{\mathrm{g}}/{\mathrm{\pi }}^2$ and ${\mathrm{R}}_0=\mathrm{R}{\mathrm{a}}_{\mathrm{g}0}/{\mathrm{\pi }}^2$ for the purposes of calculation setup and presentation. As pointed out earlier, this paper will not consider case 1 and case 5 which corresponds to the Neumann thermal boundary condition. For the Neumann thermal boundary condition, the formulation for Equation (11) needs to be changed to the appropriate function; however, that is outside the scope of the current study. For case 2 when A = 0.5 (i.e., perfectly conducting and impermeable/rigid, cooled upper and lower walls), Equation (21c) reveals that the Raleigh number is infinitely large, in other words, $\mathrm{R}\to \infty$. The results presented by [17] for a rigid lower wall and free/permeable upper wall for the same thermal boundary conditions as case 2 yields a critical internal Rayleigh number equal to 470. For the current study, case 2 shows that with internal heat generation, the presence of two rigid walls yields an infinitely large critical Rayleigh number. In practical applications, the implication of the result is that the configuration illustrated by case 2, in Table 1, indicates that the conduction solution is stable and convection will not occur.

Figure 2 shows the critical Rayleigh number versus the vibration frequency for selected values of the scaled Vadasz number for case 3 when A = −0.5 (i.e., perfectly conducting and impermeable, cooled upper walls and heated lower walls). In Figure 2, each of the curves end abruptly at specific frequency asymptotes, beyond which there are no more real solutions. These points also denote the end of the synchronous frequencies and the start of the subharmonic frequencies. The current formulation does not recover solutions corresponding to subharmonic frequencies. The work presented in [13,14] supports the point presented regarding the transition point from synchronous to subharmonic solutions. The impact of increasing the Vadasz number causes the onset of the subharmonic solutions at lower frequency values as indicated in Figure 2. At very low scaled Vadasz numbers, circa $\mathrm{\gamma }=0.01$, the onset of the subharmonic frequencies occur at very large values. In case 3, when there is no vibration, the critical Rayleigh number for internal heat generation to rank N = 1 is ${\mathrm{R}}_0=3.38{\mathrm{\pi }}^2$. This value of Rayleigh number is significantly lower than the modified internal Rayleigh number of 470 predicted by [17]. The result generated by case 3 seems plausible as the porous layer is heated from below with internal heat generation, which probably would allow for the onset of convection at a lower Raleigh number.

Figure 3 shows that on the bottom end of the vibration frequency range viz. $\mathrm{\omega }=3000$, increasing or decreasing the scaled Vadasz number has no appreciable impact on the critical Rayleigh number. However, increasing the vibration frequency to $\mathrm{\omega }=6000$ shows that increasing the scaled Vadasz number stabilizes the convection. Both Figure 2 and Figure 3 also show another important point, in other words, as the scaled Vadasz number is increased to very large values, the impact of vibration on convection becomes negligible. This observation can be recovered from the momentum equation when the time derivative disappears at very large Vadasz numbers.

Figure 4 shows the results for case 4 when A = 1.5 (i.e., perfectly conducting and impermeable, heated upper walls and cooled lower walls). In this case the results show that the existence of vibration in a system has a destabilizing effect on the convection. When the vibration frequency approaches zero in this case, the critical Rayleigh number becomes infinitely large. Even with internal heat generation, it appears that the presence of rigid/impermeable upper and lower walls results in an infinitely large Rayleigh number when the vibration frequency is zero. The results agree fully with the findings in [31,32]. In practical applications the results show that convection could occur when the vibration frequency increases to values larger than zero, in systems where the upper and lower walls are rigid/impermeable. The results in Figure 4 show that increasing the Vadasz number has a stabilizing effect on the convection. At very large Vadasz numbers, the effect of vibration becomes negligible when the contribution of the time derivative disappears in the momentum equation.

Figure 5 shows the results for the case of a porous layer heated from below and cooled from above (Benard convection) and subjected to g-jitter [13] so that the results in Figure 2, Figure 3 and Figure 4 may be compared to observe the effects of internal heat generation and the temperature boundary conditions. The impact of internal heat generation and the imposed temperature boundary conditions is to translate the curves in Figure 2 to the right when compared to the data presented in Figure 2. It is observed in all of the cases presented that the impact of increasing the Vadasz number is to have a stabilizing effect on the convection.

5. Conclusions

The results for internal heat generation with specific Dirichlet thermal boundary conditions as indicated in Table 1 are investigated when vibration is present. In this paper, case 2 to case 4 are considered. The Neumann thermal boundary conditions denoted by case 1 and 5 will be analysed separately.

When the upper and lower walls are cooled (case 2), it is demonstrated that the critical Rayleigh number is infinitely large for upper and lower rigid/impermeable walls. The work of [17] shows that for a rigid lower wall and a free/permeable upper wall, for the same thermal boundary conditions, the internal critical Rayleigh number is 470. The results for case 2 show that the basic state solution is stable and convection does not occur.

For the case of heating from below and cooling from above (case 3) it is observed that for very large Vadasz numbers the impact of vibration on convection becomes negligible. The results for case 3 are also compared to [17] and also referenced to the case of Benard convection with vibration and no internal heat generation. For case 3 it is shown that increasing the Vadasz number has a stabilizing effect on the convection. When the vibration is zero, the critical Rayleigh number approaches ${\mathrm{R}}_0=3.38{\mathrm{\pi }}^2$. This value is significantly lower than the critical internal Rayleigh number of 470 cited in [17] for cooled/free upper wall and cooled/rigid lower wall. In this case it is clear that heated lower walls combined with internal heat generated significantly lowers the convection threshold. The Vadasz number is also shown to impact the transition point from synchronous to subharmonic convection, and the results indicate that increasing the Vadasz number causes this transition point to occur sooner for case 3 and for the case of Benard convection with vibration, without internal heat generation. Upon comparing case 3 to the Benard case with vibration, there is a delay in the onset on subharmonic frequencies as shown by the translation of the curves.

When the lower walls are cooled and the upper walls are heated (case 4), it is observed that the vibration has a destabilizing effect on the convection. This result is in full agreement with previous observations [31,32].

Finally it is also demonstrated that increasing the Vadasz number stabilizes the convection in a porous layer with internal heat generation and vibration.