Nonlinear throughflow and internal heating effects on vibrating porous medium

==Abstract==

The effect of vertical throughflow and internal heating effects on fluid saturated porous medium under gravity modulation is investigated. The amplitude of modulation is considered to be very small and the disturbances are expanded in terms of power series of amplitude of convection. A weakly nonlinear stability analysis is proposed to study stationary convection. The Nusselt number is obtained numerically to present the results of heat transfer while using Ginzburg–Landau equation. The vertical throughflow has dual effect either to destabilize or to stabilize the system for downward or upward directions. The effect of internal heat source (<math display="inline">R_i>0</math>) enhances or sink (<math display="inline">R_i<0</math>) diminishes heat transfer in the system. The amplitude and frequency of modulation have the effects of increasing or diminishing heat transport. For linear model Venezian approach suggested that throughflow and internal heating have both destabilizing and stabilizing effects for suitable ranges of <math display="inline">\Omega </math>. Further, the study establishes that heat transport can be controlled effectively by a mechanism that is external to the system throughflow and gravity modulation.

==Keywords==

Throughflow; Gravity modulation; Linear and weakly nonlinear model; Internal heating

==Nomenclature==

===Latin Symbols===

<math display="inline">A</math>- amplitude of convection

<math display="inline">\delta </math>- amplitude of gravity modulation

''d''- depth of the fluid layer

<math display="inline">\overset{\rightarrow}{g}</math>- acceleration due to gravity

<math display="inline">k_c</math>- critical wave number

''K''- permeability of porous medium

''Nu''- Nusselt number

''Pe''- Péclet number <math display="inline">Pe=\frac{w_0d}{{\kappa }_T}</math>

''p''- reduced pressure

''Ra''- thermal Rayleigh number, <math display="inline">Ra=\frac{{\beta }_Tg\Delta TdK}{\nu {\kappa }_T}</math>

<math display="inline">R_i</math>- internal Rayleigh number, <math display="inline">R_i=\frac{{Qd}^2}{{\kappa }_T}</math>

<math display="inline">R_0</math>- critical Rayleigh number

''T''- temperature

<math display="inline">{Pr}_D</math>- Prandtl Darcy number <math display="inline">{Pr}_D=\frac{\phi \nu d^2}{K{\kappa }_T}</math>

<math display="inline">\Delta T</math>- temperature difference across the porous layer

''t''- time

<math display="inline">\left(x\mbox{,}z\right)</math>- horizontal and vertical co-ordinates

===Greek symbols===

<math display="inline">{\beta }_T</math>- coefficient of thermal expansion

<math display="inline">{\delta }^2</math>- square of horizontal wave number <math display="inline">{\delta }^2=k_c^2+{\pi }^2</math>

<math display="inline">\in </math>- perturbation parameter

<math display="inline">{\kappa }_T</math>- effective thermal diffusivity

<math display="inline">\gamma </math>- heat capacity ratio <math display="inline">\gamma =\frac{{\left(\rho c\right)}_m}{{\left(\rho c\right)}_f}</math>

<math display="inline">\Omega </math>- frequency of modulation

<math display="inline">\mu </math>- dynamic viscosity of the fluid

<math display="inline">\phi </math>- porosity

<math display="inline">\nu </math>- kinematic viscosity, <math display="inline">\left(\frac{\mu }{{\rho }_0}\right)</math>

<math display="inline">\rho </math>- fluid density

<math display="inline">\psi </math>- stream function

<math display="inline">\tau </math>- slow time <math display="inline">\tau ={\in }^2t</math>

<math display="inline">\theta </math>- phase angle

<math display="inline">T^{{'}}</math>- perturbed temperature

===Other symbols===

<math display="inline">{\nabla }^2</math>- <math display="inline">\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}+</math><math>\frac{{\partial }^2}{\partial z^2}</math>

===Subscripts===

''b''- basic state

''c''- critical

0- reference value

===Superscripts===

<math display="inline">{'}</math>- perturbed quantity

<math display="inline">{_\ast}</math>- dimensionless quantity

==1. Introduction==

Natural convection (buoyancy driven convection, in which gravitational force plays a main role) in fluid saturated porous medium is of interest due to its importance in various practical applications such as oil recovery process in petroleum industry, geothermal energy extraction, and insulation of reactor vessels. The above mentioned applications along with porous media occur in many natural situations. In these applications control of convective instability plays an important role. It is fact that, one of the effective mechanisms which control convective instability is to maintain a nonlinear temperature gradient. Such a gradient may be achieved by the following: time dependent heating or cooling at the boundaries, suitable thermal and rotation modulation, vibrating the porous medium periodically, volumetric distribution of internal heat sources and radiative heat transfer.

A modified complex body force (gravity modulation) is important when the system is under vertical vibrations. In this case the density gradient is subjected to vibrations, and the resulting buoyancy forces which are produced by the interaction of the density gradient with gravitational field have a complex spatiotemporal structures. The time dependent gravity field is of interest in space laboratory experiments, in areas of crystal growth and large-scale convection of atmosphere other applications. Many theoretical and experimental studies dealing with materials processing or physics of fluids under the microgravity conditions aboard an orbiting spacecraft have been carried out by Nelson [[#b0005|[1]]]. According to Wadih et al. [[#b0010|[2]]] and [[#b0015|[3]]], the vibrations can either substantially enhance or retard heat transfer and thus drastically affect the convection. The effect of modulated gravity on a convectively stable configuration can significantly influence the stability of a system by enhancing or decreasing its susceptibility to convection. Gershuni et al. [[#b0020|[4]]] and Gresho and Sani [[#b0025|[5]]] were the first to study the gravity modulation on the stability of a heated fluid layer. Their results show that the stability of the layer being heated from below is enhanced by g-jitter and being heated from below is enhanced by g-jitter. The instability of a viscoelastic fluid layer heated from below in a modulated gravitational field is studied numerically by Yang [[#b0030|[6]]]. It is found that, modulation has a destabilization effect at low frequencies and a slight stabilization effect at high frequencies, which increases with increasing the amplitude of modulation. It is also pointed that, for viscoelastic fluid acted on by a modulated gravity, modulation has the same effects at both very low and very high frequencies, as those of Newtonian fluids. Using Venezian model the onset of convection in both fluid and porous layers is investigated by Malashetty and Padmavathi [[#b0035|[7]]]. It is found that the low frequency g-jitter can have a significant effect on the stability of the problem. The Darcy limit and viscous flow limit are obtained as degenerate cases of the Brinkman model. The boundary-layer flow induced by a constant temperature vertical surface embedded in a vibrating porous medium is studied by Rees and Pop [[#b0040|[8]]]. The amplitude of g-jitter is assumed to be small compared with the mean acceleration. It is found that, the effect of g-jitter is eventually confined to a thin layer embedded within the main boundary layer, but it becomes weak at increasing distances from the leading edge. A weakly nonlinear analysis in a porous medium under gravity modulation is investigated by Govender [[#b0045|[9]]]. It is found that, the vibration frequency causes the convection amplitude to approach zero which means increasing the vibration frequency stabilizes the system. Kuznetsov [[#b0050|[10]]] investigated bioconvection in a shallow horizontal fluid-saturated porous layer that contains a suspension of oxytactic bacteria, such as Bacillus subtilis. It is found that, linear stability analysis indicates g-jitter has a stabilizing effect on the suspension. Some of the relevant and documented works on gravity modulation are [[#b0055|[11]]], [[#b0060|[12]]], [[#b0065|[13]]], [[#b0070|[14]]] and [[#b0075|[15]]].

Several geophysical and technological applications involve non-isothermal flow of fluids through porous media, called throughflow where the basic flows not quiescently. In situ processing of energy resources such as coal, geothermal energy, oil shale and many practical problems often involves the throughflow in the porous medium. The importance of buoyancy-driven convection in such problem may become significant when precise processing is required. Moreover, the throughflow effect in such circumstances may be of interest because of the possibility of controlling the convective instability by adjusting the throughflow in addition to the gravitational force due to gravity. Throughflow alters the basic state temperature from linear to nonlinear with layer height, which in turn affects the stability of the system significantly. It is found that the throughflow is stabilizing if the bounding surfaces of the porous layer are of the same type but destabilizing in one particular direction if the boundaries are not of the same type [[#b0080|[16]]], [[#b0085|[17]]], [[#b0090|[18]]] and [[#b0095|[19]]]. However, it has been shown that the throughflow is destabilizing even if the boundaries are of the same type if the porous layer is heated internally [[#b0100|[20]]] and also in the presence of an additional diffusing component. The effect of throughflow on the onset of convection in a horizontal porous layer for a Newtonian fluid has been studied by Nield [[#b0105|[21]]]. The effect of throughflow and variable gravity field on thermal convection in a porous layer is investigated by Suma et al. [[#b0110|[22]]]. The Darcy–Forchheimer model was considered for insulating boundaries. It is found that, stabilizing or destabilizing the system depends on types of boundary conditions, direction of throughflow and variable gravity field. Nonlinear throughflow and g-jitter effects on porous medium for stationary and oscillatory mode of convection are investigated by Kiran [[#b0115|[23]]] and [[#b0120|[24]]]. Weakly nonlinear model was developed in order to study heat and mass transfer across the porous layer, while deriving the non-autonomous Ginzburg–Landau equation. The similar problem for temperature modulation is studied by Kiran [[#b0125|[25]]]. Where three types of temperature profiles considered and found that, when the boundaries are at out of phase modulation the heat or mass transfer is maximum. The stationary mode nonlinear Darcy convection is investigated. Throughflow and internal heating effect on onset of convection in a porous medium is investigated by Vanishree [[#b0130|[26]]]. The effect of throughflow and internal heat generation on onset of convection is discussed. It is also found that, the onset of convection may stabilize or destabilize depending on the direction of the flow and internal heat generation.

In experiments like injection of the energy in a system is controlled by the flow rate of fluid (since the granular materials are dissipative injection of the energy is necessary to preserve steady states). The system is in the fixed phase at low flow rate, where particle rests at the bottom. When the flow rate exceeds the critical value, the particles start moving, this state is called the fluidized phase. These fluidized beds may be employed and are created by throughflow through solid particles. Throughflow is important for industrial applications in that it yields a potential method to control the onset behavior of convection, since many practical problems involve more than a Boussinesq approximation in the system of equations. The aim of throughflow is to develop a nonlinear energy stability analysis which yields unconditional stability bounds which are closed to linear theory. It is quite important for understanding these nonlinear effects in nonlinear thermal instability to finite amplitude convection.

The study of heat transfer with internal heat generation arises at consideration of heat dissipation in porous layers exposed to radioactive radiation at storage of nuclear waste materials or transmitting high energy particle beams. There are several situations where porous material offers its own heat source and set up convective flow through the local heat generation. This kind of problem may occur through radioactive decay, a relatively weak exothermic reaction which can take place within the porous material. To be specific, it is the main source of energy for celestial bodies caused by nuclear fusion and decaying of radioactive materials, which keeps them warm and active. Therefore, internal heating of is important in several applications that include reactor safety analyses, geophysics, fire and combustion studies, metal waste form development for spent nuclear fuel, and nuclear reactions, nuclear heat cores, nuclear energy, nuclear waste disposals, oil extractions and crystal growth. The heat transfer in a window for transmission of high-energy beams with a metal foam porous layer is studied by Hetsroni et al. [[#b0135|[27]]]. The image processing of thermal maps on the surface of open-cell porous layer revealed that, the boundary condition with constant heat flux became more relevant for numerical model than those with constant temperature at high values of Reynolds numbers based on the permeability. It is also found that, heat transfer at natural convection in the strip of metal foam was increased drastically (up to 1820 times for metal foam of 20 ppi) relative to the flat plate of the same overall dimensions. Bhadauria et al. [[#b0140|[28]]] performed a weakly nonlinear thermal instability in a porous medium in the presence of internal heating and gravity modulation. An amplitude of convection is derived using Landau equation and found gravity modulation used, to alter the heat transfer, whereas internal heating is to enhance the heat transfer. The studies related to internal heat generation were recently investigated by Bhadauria et al. [[#b0145|[29]]] and Saravanan [[#b0150|[30]]].

It is observed that, in the literature many models being implemented with throughflow for linear theory. But, for nonlinear models with modulation it is scary. It is due to Kiran et al. [[#b0115|[23]]], [[#b0120|[24]]] and [[#b0125|[25]]] that, throughflow effects are considered under gravity or temperature modulation. These are the first papers which contribute the nonlinear throughflow effects under modulation using stationary and oscillatory mode of convections. The conduction state temperature is nonlinear in terms of Péclet number which measures the strength of throughflow and it affects the system through coupled equations of momentum and energy. It is found that, this nonlinear throughflow has duel effect on heat transfer in the system. Also three types of temperature modulation are discussed and presented a significant results for out of phase modulation. A linear stability analysis is performed by Vanishree [[#b0130|[26]]] to study the effect of throughflow and internal heat generation on the preferred mode of stationary thermal convection in a variable viscosity liquid saturating an anisotropic porous medium. The results are presented corresponding to the onset of convection. But, unfortunately the nonlinear part is missing therefore heat transfer analysis. With this concussion the present article is motivated to continue the previous works of Kiran et al. [[#b0125|[25]]] and [[#b0195|[39]]] and Vanishree [[#b0130|[26]]] for nonlinear model. Hence the present paper contributes nonlinear throughflow effects under gravity modulation in the presence of internal heat source and sink, while employing weakly nonlinear model.

==2. Mathematical formulation==

Consider an infinitely extended horizontal porous medium saturated by Newtonian fluid, confined between two free-free boundaries at <math display="inline">z=0</math> and <math display="inline">z=d</math>, and heated from below. The temperature difference across the porous medium is kept at <math display="inline">\Delta T</math>. The Cartesian frame of reference is chosen, origin in the lower boundary and the ''z'' axis in vertically upward direction. The schematic diagram is given in the  [[#f0005|Fig. 1]]. Throughflow has been considered in vertical, both upward and downward directions. Further Darcy law and the Oberbeck–Boussinesq approximation are considered, under these assumptions, the equations which describe the flow model are given by Bhadauria and Kiran [[#b0155|[31]]]:

<span id='e0195'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\nabla .\overset{\rightarrow}{q}=0\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

<span id='e0200'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{{\rho }_0}{\phi }\frac{\partial \overset{\rightarrow}{q}}{\partial t}=</math><math>-\nabla p+\rho \overset{\rightarrow}{g}-\frac{\mu }{K}\overset{\rightarrow}{q}\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

<span id='e0205'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\gamma \frac{\partial T}{\partial t}+(\overset{\rightarrow}{q}\cdot \nabla )T=</math><math>{\kappa }_T{\nabla }^2T+Q(T-T_0)\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

<span id='e0210'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\rho ={\rho }_0\left[1-{\beta }_T\left(T-T_0\right)\right]\mbox{.}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

The physical variable has their usual meanings and is given in Nomenclature. For simplicity thermal expansion coefficient <math display="inline">\gamma </math> is taken to be unity in this paper. The externally imposed thermal boundary conditions and time periodic gravity field considered in this paper are:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{rl}
T & =T_0+\frac{\Delta T}{2}\mbox{,}\quad \mbox{at}\quad z=0\\
 & =T_0\mbox{,}\quad \mbox{at}\quad z=d
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\overset{\rightarrow}{g}=g_0(1+\delta {\in }^2cos(\Omega t))\hat{k}\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

where <math display="inline">\delta </math> and <math display="inline">\Omega </math> represent the amplitude and frequency of gravity modulation, respectively. The basic state is assumed to be quiescent and the quantities in this state are given by:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>q_b=(0\mbox{,}0\mbox{,}w_0)\mbox{,}\quad \rho ={\rho }_b(z\mbox{,}t)\mbox{,}\quad p=</math><math>p_b(z\mbox{,}t)\mbox{,}\quad T=T_b(z)\mbox{.}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial p_b}{\partial z}=\frac{\mu }{K}w_0-</math><math>{\rho }_bg\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

<span id='e0015'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>w_0\frac{\partial T_b}{\partial z}={\kappa }_T\frac{{\partial }^2T_b}{\partial z^2}+</math><math>Q(T_b-T_0)\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\rho }_b={\rho }_0\left[1-{\beta }_T\left(T_b-T_0\right)\right]\mbox{.}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

The basic state temperature <math display="inline">T_b</math> defined numerically by solving the Eq. [[#e0015|(9)]]. The finite amplitude perturbations on the basic state are superposed in the form,

<span id='e0025'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\overset{\rightarrow}{q}={\overset{\rightarrow}{q}}_b+</math><math>{\overset{\rightarrow}{q}}^{{'}}\mbox{,}\quad \rho =</math><math>{\rho }_b+{\rho }^{{'}}\mbox{,}\quad p=p_b+p^{{'}}\mbox{,}\quad T=</math><math>T_b+T^{{'}}\mbox{.}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

For two dimensional flow model the following stream function introduced <math display="inline">\psi </math> as <math display="inline">u^{{'}}=\frac{\partial \psi }{\partial z}\mbox{,}w^{{'}}=</math><math>-\frac{\partial \psi }{\partial x}</math>. The non-dimensional physical variables re-scaled as <math display="inline">(x\mbox{,}y\mbox{,}z)=d(x^{{_\ast}}\mbox{,}y^{{_\ast}}\mbox{,}z^{{_\ast}})\mbox{,}(u\mbox{,}w)=</math><math>\frac{{\kappa }_T}{d}(u^{{_\ast}}\mbox{,}w^{{_\ast}})\mbox{,}p=</math><math>\frac{\mu {\kappa }_T}{K}p^{{_\ast}}\mbox{,}\Omega =</math><math>\frac{{\kappa }_T}{d^2}{\Omega }^{{_\ast}}</math><math display="inline">\quad t=\frac{d^2}{{\kappa }_T}t^{{_\ast}}\mbox{,}\quad \psi =</math><math>{\kappa }_T{\psi }^{{_\ast}}\mbox{,}\quad T^{{'}}=</math><math>\Delta T\quad T^{{_\ast}}</math>. Substituting the Eq. [[#e0025|(11)]] in Eqs.  [[#e0195|(1)]], [[#e0200|(2)]], [[#e0205|(3)]] and [[#e0210|(4)]], and using the basic state results, while eliminating the pressure term, the following dimensionless governing equations arrived (dropping the asterisk):

<span id='e0030'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{1}{{Pr}_D}\frac{\partial }{\partial t}({\nabla }^2\psi )=</math><math>-{\nabla }^2\psi -Ra\quad g_m\frac{\partial T}{\partial x}\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

<span id='e0035'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>-\frac{{dT}_b}{dz}\frac{\partial \psi }{\partial x}-</math><math>\left({\nabla }^2+R_i-Pe\frac{\partial }{\partial z}\right)T=</math><math>-\frac{\partial T}{\partial t}+\frac{\partial (\psi \mbox{,}T)}{\partial (x\mbox{,}z)}\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

where <math display="inline">g_m=(1+{\in }^2\delta cos(\Omega t))</math>. The dimensionless parameters which appear in Eq.  [[#e0030|(12)]] and [[#e0035|(13)]] are given in Nomenclature. The Eq. [[#e0035|(13)]] shows the basic state solution and influences the stability problem through the factor <math display="inline">\frac{{dT}_b}{dz}</math>, which is given by Eq. [[#e0040|(14)]]. Assuming small variation in time, and re-scaling it as <math display="inline">\tau ={\in }^2t</math>, the stationary mode of convection is studied. The basic state temperature appears in Eq. [[#e0035|(13)]] is derived from the Eq. [[#e0015|(9)]] numerically and it is given by:

<span id='e0040'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Pe\frac{{dT}_b}{dz}=\frac{d^2T_b}{{dz}^2}+R_iT_b\mbox{.}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

The nonlinear system of Eqs.  [[#e0030|(12)]] and [[#e0035|(13)]] may be written in the following matrix form:

<span id='e0045'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left[\begin{array}{cc}
{\nabla }^2 & g_mRa\frac{\partial }{\partial x}\\
-\frac{\partial T_b}{\partial z}\frac{\partial }{\partial x} & -\left({\nabla }^2+R_i-Pe\frac{\partial }{\partial z}\right)
\end{array}\right]\left[\begin{array}{c}
\psi \\
T
\end{array}\right]=\left[\begin{array}{c}
-\frac{{\in }^2}{{Pr}_D}\frac{\partial }{\partial \tau }({\nabla }^2\psi )\\
-{\in }^2\frac{\partial T}{\partial \tau }+\frac{\partial (\psi \mbox{,}T)}{\partial (x\mbox{,}z)}
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

<span id='f0005'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_García_631601349-1-s2.0-S1110016816000181-gr1.jpg|center|369px|A schematic diagram of the problem.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figure 1.

A schematic diagram of the problem.

</span>
|}

The solution of the system Eq. [[#e0045|(15)]], is evaluated by considering impermeable boundary conditions as given (Bhadauria and Kiran [[#b0155|[31]]]):

<span id='e0225'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{rl}
 & \psi =0\quad \mbox{and}\quad T=0\quad \mbox{at}\quad z=0\mbox{,}\\
 & \psi =0\quad \mbox{and}\quad T=0\quad \mbox{at}\quad z=1\mbox{.}
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

==3. Heat transport for stationary instability==

The solution of the above system Eq. [[#e0045|(15)]] is derived using an asymptotic expansions (given by Malkus and Veronis [[#b0160|[32]]]):

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{rl}
 & Ra=R_0+{\in }^2R_2+{\in }^4R_4+\cdots \mbox{,}\\
 & \psi =\in {\psi }_1+{\in }^2{\psi }_2+{\in }^3{\psi }_3+\cdots \mbox{,}\\
 & T=\in T_1+{\in }^2T_2+{\in }^3T_3+\cdots \mbox{,}
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

where <math display="inline">R_0</math> is the critical value of the Rayleigh number at which the onset of convection takes place in the absence of gravity modulation. The above expansion was earlier used by Malkus and Veronis [[#b0160|[32]]] and Venezian [[#b0165|[33]]]. Now the solution of the system evaluated for different orders of <math display="inline">\in </math>.

===3.1. Lowest-order system===

This lowest order system is similar to linear model where the nonlinear effects do not occur. At this order the system takes the following form:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left[\begin{array}{cc}
{\nabla }^2 & R_0\frac{\partial }{\partial x}\\
-\frac{{dT}_b}{dz}\frac{\partial }{\partial x} & -\left({\nabla }^2+R_i-Pe\frac{\partial }{\partial z}\right)
\end{array}\right]\left[\begin{array}{c}
{\psi }_1\\
T_1
\end{array}\right]=\left[\begin{array}{c}
0\\
0
\end{array}\right]\mbox{.}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

The solution of the lowest order system is subject to the boundary conditions Eq. [[#e0225|(16)]] is:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\psi }_1=A(\tau )sin(k_cx)sin(\pi z)\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>T_1=\frac{2k_cI_1}{\left({\delta }^2-R_i\right)}A(\tau )cos(k_cx)sin(\pi z)\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

where <math display="inline">{\delta }^2=k_c^2+{\pi }^2</math> and <math display="inline">I_1={\int }_0^1\frac{{dT}_b}{dz}{sin}^2(\pi z)dz</math>. The critical value of the Rayleigh number for stationary convection is given by

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>R_0=\frac{{\delta }^2(R_i-{\delta }^2)}{2k_c^2I_1}\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|}

and the corresponding wavenumber for the onset of convection is obtained while minimizing <math display="inline">R_0</math> with respect to <math display="inline">k_c^2</math>.

===3.2. Second-order system===

In this order the nonlinear effects enter through the Jacobian terms where the interaction of the fluid motion and thermal diffusivity take place. At this order the system take the following form:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left[\begin{array}{cc}
{\nabla }^2 & R_0\frac{\partial }{\partial x}\\
-\frac{\partial T_b}{\partial z}\frac{\partial }{\partial x} & -\left({\nabla }^2+R_i-Pe\frac{\partial }{\partial z}\right)
\end{array}\right]\left[\begin{array}{c}
{\psi }_2\\
T_2
\end{array}\right]=\left[\begin{array}{c}
R_{21}\\
R_{22}
\end{array}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|}

where

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>R_{21}=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>R_{22}=\frac{\partial {\psi }_1}{\partial x}\frac{\partial T_1}{\partial z}-</math><math>\frac{\partial {\psi }_1}{\partial z}\frac{\partial T_1}{\partial x}\mbox{.}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|}

The second order solutions subjected to the boundary conditions Eq. [[#e0225|(16)]] is given by:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\psi }_2=0\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>T_2=\frac{2k_c^2\pi I_1}{\left(4{\pi }^2-R_i)({\delta }^2-R_i\right)}A^2(\tau )sin(2\pi z)\mbox{.}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

The horizontally averaged Nusselt number ''Nu'', for the stationary mode of convection is determined numerically by using the definition

<span id='e0100'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Nu(\tau )=1+\frac{{\left[\frac{k_c}{2\pi }{\int }_0^{\frac{2\pi }{k_c}}\left(\frac{\partial T_2}{\partial z}\right)\quad dx\right]}_{z=0}}{{\left[\frac{k_c}{2\pi }{\int }_0^{\frac{2\pi }{k_c}}\left(\frac{{dT}_b}{dz}\right)\quad dx\right]}_{z=0}}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

In the above Eq. [[#e0100|(27)]] 1 indicates, heat transport in conduction form and second term represents the heat transfer in convection form.

===3.3. Third-order system===

At this order the system takes the following form:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\left[\begin{array}{cc}
{\nabla }^2 & R_0\frac{\partial }{\partial x}\\
-\frac{\partial T_b}{\partial z}\frac{\partial }{\partial x} & -\left({\nabla }^2+R_i-Pe\frac{\partial }{\partial z}\right)
\end{array}\right]\left[\begin{array}{c}
{\psi }_3\\
T_3
\end{array}\right]=\left[\begin{array}{c}
R_{31}\\
R_{32}
\end{array}\right]\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

where

<span id='e0110'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>R_{31}=-\frac{1}{{Pr}_D}\frac{\partial }{\partial \tau }({\nabla }^2{\psi }_1)-</math><math>R_0\delta cos(\Omega \tau )\frac{\partial T_1}{\partial x}-</math><math>R_2\frac{\partial T_1}{\partial x}\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

<span id='e0115'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>R_{32}=\frac{\partial {\psi }_1}{\partial x}\frac{\partial T_2}{\partial z}-</math><math>\frac{\partial T_1}{\partial \tau }\mbox{.}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

Substituting the expressions of <math display="inline">{\psi }_1\mbox{,}T_1</math> and <math display="inline">T_2</math> into Eqs. [[#e0110|(29)]] and [[#e0115|(30)]], the expressions of <math display="inline">R_{31}</math> and <math display="inline">R_{32}</math> easily simplified. But, finding the solution at this order is very difficult due to complexity. The solvability condition being used for the existence of third order solution, and consequently the Ginzburg–Landau amplitude equation is derived for stationary mode of convection with time-periodic coefficients in the following form:

<span id='e0120'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>A_1\frac{dA(\tau )}{d\tau }+A_2A(\tau )-A_3A^3(\tau )=</math><math>0\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}

where <math display="inline">A_1=\frac{{\delta }^2}{{Pr}_D}-\frac{2R_0k_c^2I_1}{{\left({\delta }^2-R_i\right)}^2}</math>,  <math display="inline">A_2=\frac{2k_c^2I_1}{\left({\delta }^2-R_i\right)}(R_2+</math><math>R_0\delta cos(\Omega \tau ))\mbox{,}A_3=\frac{4R_0{\pi }^2k_c^4I_1}{2{\left({\delta }^2-R_i\right)}^2(4{\pi }^2-R_i)}</math>. The Ginzburg–Landau equation given by Eq. [[#e0120|(31)]] is a Bernoulli equation and obtaining its analytical solution is difficult due to its non-autonomous nature. So that it has been solved numerically using the in-built function NDSolve of mathematica subjected to the initial condition <math display="inline">A(0)=a_0</math>, where <math display="inline">a_0</math> is the chosen initial amplitude of convection. Without loss of generality <math display="inline">R_2=R_0</math>, is assumed in the calculations and this is done to keep the parameters to a minimum. For unmodulated case, the analytical solution of the above Eq. [[#e0120|(31)]] takes the form:

<span id='e0125'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>A(\tau )=\frac{1}{\sqrt{\left(\frac{A_3}{A_2}+C_1\mbox{Exp}\left[-\frac{2A_2}{A_1}\tau \right]\right)}}\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

where <math display="inline">A_1\mbox{,}A_3</math> same as in Eq. [[#e0120|(31)]], <math display="inline">A_2=\frac{2k_c^2I_1R_0}{\left({\delta }^2-R_i\right)}</math> and <math display="inline">C_1</math> which appears in Eq. [[#e0125|(32)]] is an integration constant, can be found by using suitable initial condition.

==4. Linear stability==

For linear stability analysis the nonlinear terms excluded from the Eqs. [[#e0030|(12)]] and [[#e0035|(13)]] thus obtain the linear model in the following form:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\nabla }^2\left(\frac{1}{{Pr}_D}\frac{\partial }{\partial t}-\right. </math><math>\left. 1\right)\psi =Ra\quad g_m\frac{\partial T}{\partial x}\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>I_B\frac{\partial \psi }{\partial x}=\left(\frac{\partial }{\partial t}-\right. </math><math>\left. {\nabla }^2+R_i+Pe\frac{\partial }{\partial z}\right)T\mbox{.}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

Eliminating ''T  '' from the above equations we get single equation for <math display="inline">\psi </math> in the form:

<span id='e0140'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\nabla }^2\left(\frac{1}{{Pr}_D}\frac{\partial }{\partial t}-\right. </math><math>\left. 1\right)\left(\frac{\partial }{\partial t}-\right. </math><math>\left. {\nabla }^2+R_i+Pe\frac{\partial }{\partial z}\right)\psi -</math><math>I_B\quad Ra\quad g_m\frac{{\partial }^2\psi }{\partial x^2}=</math><math>0\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}

where the gravity modulation term takes the form <math display="inline">g_m=1+\in cos(\Omega t)</math> for linear case.

===4.1. Method of solution===

We obtain the perturbation technique to obtain eigen functions <math display="inline">\psi </math> and eigen values Ra by introducing a small perturbation parameter <math display="inline">\in </math> (Venezian [[#b0165|[33]]]).

<span id='e0145'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\begin{array}{l}
\psi ={\psi }_0+\in {\psi }_1+{\in }^2{\psi }_2+\cdots \\
Ra={Ra}_0+{\in }^2{Ra}_2+\cdots 
\end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

Substituting the equation Eq. [[#e0145|(36)]] into equation Eq. [[#e0140|(35)]] the system will be solved for every order of <math display="inline">\in </math>.

===4.2. Zeroth order system===

The marginally stable solutions will be considered of the following form

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\psi }_0=sin(ax)sin(\pi z)\mbox{.}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

Since <math display="inline">\psi </math> is independent of ''y'' therefore ''a'' is the wave number in ''x'' direction, the corresponding eigen values are given by

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{Ra}_0=\frac{{\delta }^2({\delta }^2+R_i)}{a^2I_B}\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

where <math display="inline">{\delta }^2={\pi }^2+a^2</math> and <math display="inline">I_B={\int }_0^1\frac{{dT}_b}{dz}{sin}^2(\pi z)dz</math>. In the absence of internal heat, the expression of <math display="inline">{Ra}_0</math> similar with the result of Kiran et al.  [[#b0120|[24]]] and [[#b0125|[25]]].

===4.3. First order system===

At this order we obtain the following relation:

<span id='e0160'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>L({\psi }_1)=I_B{Ra}_0\in cos(\Omega t)\frac{{\partial }^2\psi }{\partial x^2}\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

where <math display="inline">L={\nabla }^2\left(\frac{1}{{Pr}_D}\frac{\partial }{\partial t}-\right. </math><math>\left. 1\right)\left(\frac{\partial }{\partial t}-\right. </math><math>\left. {\nabla }^2+R_i+Pe\frac{\partial }{\partial z}\right)-</math><math>\frac{{\partial }^2}{\partial x^2}</math>. The solution of the homogeneous equation corresponding to the equation Eq. [[#e0160|(39)]], one must have a term which is orthogonal to <math display="inline">sin(\alpha x)sin(\pi z)</math>. Note that we already apply the theory of minimal Rayleigh number therefore right hand side of equation Eq. [[#e0160|(39)]] is complete orthogonal function to <math display="inline">sin(\alpha x)sin(\pi z)</math>. Using L operator it is found that

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>L[sin(\pi z){exp}^{\left(iax+i\pi z-i\Omega \tau \right)}]=</math><math>L(\Omega )sin(\pi z){exp}^{\left(iax+i\pi z-i\Omega \tau \right)}\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>L[\Omega ]=-{\delta }^2\left(\frac{1}{{Pr}_D}(-i\Omega )-\right. </math><math>\left. 1\right)\left(-i\Omega +{\delta }^2+R_i-Pei\pi \right)+</math><math>I_B{Ra}_0a^2\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}

The solution for <math display="inline">{\psi }_1</math> will be derived by inverting the operator ''L'' in the following form

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\psi }_1=-I_B{Ra}_0a^2{\psi }_0cos(\Omega t)Re\left[\frac{L^{{_\ast}}(\Omega )}{\vert L(\Omega )\vert ^2}\right]\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
|}

where <math display="inline">L^{{_\ast}}(\Omega )</math> is the complex conjugate of <math display="inline">L(\Omega )</math>.

===4.4. Second order system===

At this order the following relation is obtained:

<span id='e0180'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>L({\psi }_2)=-I_B{Ra}_2a^2{\psi }_0+I_B{Ra}_0a^2\in cos(\Omega t){\psi }_1\mbox{.}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}

We will not solve the Eq. [[#e0180|(43)]], but using this equation we determine <math display="inline">{Ra}_2</math>. For the existence of a solution of the Eq. [[#e0180|(43)]], it is necessary that, the steady part of its right hand side must orthogonal to <math display="inline">sin(\alpha x)sin(\pi z)</math>. This gives,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{\int }_0^1{\int }_0^{\frac{\pi }{\alpha }}\left[\lbrace -\right. </math><math>\left. I_B{Ra}_2a^2{\psi }_0+{\left(I_B{Ra}_0a^2\in \right)}^2\rbrace \overline{{cos}^2(\Omega t)}\right]\frac{Re[L^{{_\ast}}(\Omega )]}{\vert L(\Omega )\vert ^2}dxdz\mbox{,}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}

where the upper bar denotes the time average,

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>{Ra}_2=\frac{-I_B{\left({Ra}_0a\in \right)}^2A}{2{\left(A^2+B^2\right)}^2}\mbox{.}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
|}

==5. Results and discussion==

The combined effect of gravity modulation and vertical throughflow is investigated in the presence of internal heat source and sink. Both linear and weakly nonlinear theories are investigated. A weakly nonlinear study investigates heat transfer analysis and linear model investigates onset convection in the medium. Firstly the result related nonlinear theory will be discussed. The gravity modulation is assumed to be of third order <math display="inline">O({\in }^2)</math>, which means only small amplitude of gravity modulation is considered. This assumption will help us in obtaining the amplitude equation of convection in simple and elegant manner and is much easier to obtain than Lorenz model. The purpose of weakly nonlinear theory is to study heat transfer, which is not possible in linear theory. External regulations (like thermal, gravity, rotation, and magnetic field modulations) of the current problem is important due to its numerous applications in thermal and engineering sciences. The objective of the problem is to consider, gravity modulation and vertical throughflow for either enhancing or diminishing heat transfer. Since the porous medium is assumed to be closely packed, the Darcy-model is considered in momentum equation. The additional important topic is internal heat source or sink which was realized by energy equation. This concept is most important where the system offers its own internal heat generation. The moderate values of <math display="inline">R_i</math> are considered due to its dominant effect on convection.

The values of <math display="inline">{Pr}_D</math> may consider (modern porous medium applications, such as mushy layer in solidification of binary alloys and fractured porous medium) around 1, and also for low porosity medium, the large values considered for <math display="inline">{Pr}_D</math>. The values of <math display="inline">\delta </math> and <math display="inline">\Omega </math> are considered to be small. For small values of amplitude and frequencies, the heat transfer can reach the maximum. The numerical results for Nu obtained from the Eq. [[#e0100|(27)]] in connection with amplitude Eq. [[#e0120|(31]]), and the results have been presented in [[#f0010|Figure 2]], [[#f0015|Figure 3]] and [[#f0020|Figure 4]]. The effect of each parameter on heat transport is shown in [[#f0010|Figure 2]], [[#f0015|Figure 3]] and [[#f0020|Figure 4]] where the plots of Nusselt number Nu versus <math display="inline">\Omega </math> are presented. It is found from the figures that the value of Nu oscillates maximum for low frequencies and further values of <math display="inline">\Omega </math> no oscillations can be observed, which means that, low frequencies of modulation enable the system to convect faster. The oscillatory nature of the figures is due to modulation only.

<span id='f0010'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_García_631601349-1-s2.0-S1110016816000181-gr2.jpg|center|px|Effect of Pe and Ri on heat transport.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figure 2.

Effect of Pe and <math display="inline">R_i</math> on heat transport.

</span>
|}

<span id='f0015'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_García_631601349-1-s2.0-S1110016816000181-gr3.jpg|center|px|Effect of PrD and δ on heat transport.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figure 3.

Effect of <math display="inline">{Pr}_D</math> and <math display="inline">\delta </math> on heat transport.

</span>
|}

<span id='f0020'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_García_631601349-1-s2.0-S1110016816000181-gr4.jpg|center|303px|Comparison between modulated and unmodulated system.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figure 4.

Comparison between modulated and unmodulated system.

</span>
|}

Let us discuss the results of gravity modulation: in [[#f0010|Fig. 2]] the results corresponding to Péclet and internal Rayleigh number present. The effect of ''Pe  '' on heat transfer is investigated for the cases of upward and downward directed flows. The upward throughflow (also known as pro-gravity) <math display="inline">\left(Pe>0\right)</math>, given in [[#f0010|Fig. 2]]a has destabilizing effect whereas downward throughflow (also known as anti-gravity) <math display="inline">\left(Pe<0\right)</math>, given in [[#f0010|Fig. 2]]b has stabilizing effect. The related results for linear theory are due to Vanishree [[#b0130|[26]]], this is because of temperature dependence of viscosity. When ''Pe'' varies in either directions the nonlinear temperature distribution significantly matters at conduction state and there is significant or drastic variation in the energy supply to the disturbance. The corresponding results are also obtained by Nield  [[#b0105|[21]]] in the case of fluid layer for small amount of throughflow. The present results are computable with the results obtained by Shivakumara et al. [[#b0095|[19]]] and Suma et al. [[#b0110|[22]]]. Shivakumara et al. [[#b0095|[19]]] pointed that, the destabilization effect may be due to the distortion of basic state temperature distribution from linear to nonlinear throughflow. The corresponding results of throughflow under modulation the reader may look at the articles  [[#b0125|[25]]] and [[#b0130|[26]]]. It is noted that, for nonlinear theories ''Pe'' has duel effect  [[#b0125|[25]]] but, for linear case <math display="inline">\vert Pe\vert </math> only stabilizing effect [[#b0130|[26]]], the reason could be guessed as the finite amplitude interaction in the fluid flow through coupled formation of momentum and energy equations.

The effect of the internal Rayleigh number <math display="inline">R_i</math>, for both heat source and heat sink cases, on heat transfer presented in [[#f0010|Fig. 2]]c for heat source and [[#f0010|Fig. 2]]d for heat sink. The reason behind internal heat is that, for many industrial and chemical experiments during their process the system may release or absorb energy in such a situation one needs understand how to control these fluctuations in the system. It is given figures that positive values of <math display="inline">R_i</math> (heat source) increase or negative values of <math display="inline">R_i</math> (heat sink) decreases the heat transfer in the system. The reader may refer [[#b0065|[13]]], [[#b0130|[26]]], [[#b0135|[27]]], [[#b0140|[28]]], [[#b0145|[29]]], [[#b0150|[30]]] and [[#b0155|[31]]] for internal heat generation. In Astrophysics provides an example of internally driven convection for uniformly heated medium. In the cores of stars heat is produced by thermonuclear reactions. In such situation the heating rate is very sensitive to temperature, and this sensitivity creates steep thermal gradients that drive powerful convection in the cores (Kippenhahn et al. [[#b0180|[36]]]). Clearly, this and many other instances (related space science) of internally driven convection contain more complications. It is observed that, for low porosity medium, large values [[#b0145|[29]]] considered for <math display="inline">{Pr}_D</math>. The Nusselt number Nu increases with <math display="inline">{Pr}_D</math> showing heat transfer increases (see [[#f0015|Fig. 3]]a). The results could be gained for lower values of Prandtl Darcy number. The related articles which presents the results corresponding to heat transfer in a porous medium under modulation may be observed in the following studies [[#b0060|[12]]], [[#b0065|[13]]] and [[#b0155|[31]]].

From the [[#f0010|Figure 2]], [[#f0015|Figure 3]] and [[#f0020|Figure 4]], it is observed that as increase <math display="inline">\Omega </math> decreases the magnitude of Nu, and so the effect of frequency of modulation on heat transport diminishes heat transfer. At high rates of frequency, the gravity modulation on thermal instability disappears altogether. The above results agree quite well with the linear theory of Venezian [[#b0165|[33]]] for temperature modulation, where the correction in the critical value of Rayleigh number due to thermal modulation becomes almost zero at high frequencies. The results for gravity modulation on weakly nonlinear studies see [[#b0035|[7]]], [[#b0055|[11]]], [[#b0115|[23]]], [[#b0120|[24]]], [[#b0125|[25]]], [[#b0130|[26]]], [[#b0135|[27]]], [[#b0140|[28]]] and [[#b0145|[29]]]. In general the frequency of modulation shortens the wavelength and minimizes the magnitude of Nu, due to this reason heat transfer reduces. It is quite intrusting to see even this frequency of modulation play critical role in controlling chaos [[#b0120|[24]]]. The effect of amplitude of modulation is presented in [[#f0015|Fig. 3]]b, and it is found that heat transfer increases as <math display="inline">\delta </math> increases. Hence amplitude of modulation is to enhance the heat transfer. The comparison between modulated and unmodulated case is presented in [[#f0020|Fig. 4]]. It is found that gravitationally modulated system flows transport less heat transfer than unmodulated systems the corresponding results obtained from the studies of [[#b0170|[34]]] and [[#b0175|[35]]]. The effect of magnetic field [[#b0185|[37]]] modulation and rotational speed [[#b0190|[38]]] modulation on thermal instability shows similar results for modulated systems. It is to be noted that, the analytical solution for unmodulated case obtained from the Eq. [[#e0125|(32)]].

Finally the nature of stream lines and corresponding isotherms presented in [[#f0025|Figure 5]] and [[#f0030|Figure 6]]. [[#f0025|Fig. 5]] shows the variation of stream lines and isotherms at different instants of times, respectively. [[#f0025|Fig. 5]] shows the magnitudes of stream lines increase as time increases. Also, initially the isotherms are flat and parallel, thus heat transport is due to conduction only. However, as time increases, isotherms form contours, showing convective regime is taking place, after reaching certain instant there is no change in the magnitude of stream lines and isotherms, thus showing the steady state. In particular [[#f0030|Fig. 6]] is drawn to see the effect of throughflow on isotherms. While strengthening the throughflow there forms a boundary layer at the bottom plate, then the significant results may observe at corresponding isotherms, where isotherms are missing at the bottom plate due to throughflow disturbances. The reason would be, when there is flow at the bottom plate there may not be stream line flow due to heavy boundary layer and hence the nature of the [[#f0030|Fig. 6]].

<span id='f0025'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_García_631601349-1-s2.0-S1110016816000181-gr5.jpg|center|px|Streamlines and Isotherms for Pe=0.3, Ri=1.0PrD=100∊=0.5.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figure 5.

Streamlines and Isotherms for Pe = 0.3, <math display="inline">R_i=1.0</math><math display="inline">{Pr}_D=100</math><math display="inline">\in =0.5</math>.

</span>
|}

<span id='f0030'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_García_631601349-1-s2.0-S1110016816000181-gr6.jpg|center|446px|Isothermas for τ=4.0Ri=1.0PrD=100∊=0.5.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figure 6.

Isothermas for <math display="inline">\tau =4.0</math><math display="inline">R_i=1.0</math><math display="inline">{Pr}_D=100</math><math display="inline">\in =0.5</math>.

</span>
|}

In order to make the article self contained the results corresponding to linear theory also presented in [[#f0035|Fig. 7]]. It is examined that, how this correction Rayleigh number <math display="inline">{Ra}_2</math> affected by the modulation frequency <math display="inline">\Omega </math>. Analytical study of gravity modulation on onset of thermal convection in a porous layer using Venezian approach is presented. The regular perturbation method based on small amplitude of modulation is employed to compute the critical value of Rayleigh number and wavenumber. The expression for the critical correction Rayleigh number <math display="inline">{Ra}_2</math> is computed as a function of the frequency of modulation <math display="inline">\Omega </math>. The sign of <math display="inline">{Ra}_2</math> characterizes the stabilizing or destabilizing effect of modulation. A positive <math display="inline">{Ra}_2</math> indicates that the modulation effect is stabilizing while a negative <math display="inline">{Ra}_2</math> indicates destabilizing effect.

<span id='f0035'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_García_631601349-1-s2.0-S1110016816000181-gr7.jpg|center|556px|Variations of Ra2 with Ω for different values of PrD, Pe and Ri.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figure 7.

Variations of <math display="inline">{Ra}_2</math> with <math display="inline">\Omega </math> for different values of <math display="inline">{Pr}_D</math>, Pe and <math display="inline">R_i</math>.

</span>
|}

In [[#f0035|Fig. 7]]a the effect of DarcyPrandtl number <math display="inline">{Pr}_D</math> on the stability of the system is displayed. It is found that the magnitude of critical correction Rayleigh number <math display="inline">{Ra}_2</math> decreases positively with increasing ''PrD'' when frequency is small. However, the trend reverses for moderate and high frequency. Therefore, increase of DarcyPrandtl number enhances the destabilizing effect of gravity modulation when frequency is small whereas it enhances the stabilizing effect of modulation for moderate and high frequency. The results of the present model may be compared with the studies of Malashetty and Padmavathi  [[#b0035|[7]]], Vanishree [[#b0130|[26]]] and Kiran et al. [[#b0195|[39]]]. [[#f0035|Fig. 7]]b presents result corresponding to internal Rayleigh number <math display="inline">R_i</math>, as <math display="inline">R_i</math> increases positively (heat source), <math display="inline">\vert {Ra}_2\vert </math> decreases positively when <math display="inline">R_i=1</math> and 2. When <math display="inline">R_i>2</math> the value of <math display="inline">\vert {Ra}_2\vert </math> increases negatively show destabilizing effect. But, for heat sink ([[#f0035|Fig. 7]]c) i.e. when <math display="inline">Ri<0</math> the magnitude of <math display="inline">{Ra}_2</math> increases showing stabilizing effect of the system. In [[#f0035|Fig. 7]]d and e, the effect of throughflow is presented. For up-directed flow <math display="inline">Pe>0</math> the value of <math display="inline">\vert {Ra}_2\vert </math> decreases positively showing destabilizing effect. While for inflow the value of <math display="inline">\vert {Ra}_2\vert </math> negative for lower values of <math display="inline">\Omega </math> and decreasing the value of <math display="inline">\vert {Ra}_2\vert </math> showing stabilizing effect. Further moderate values of <math display="inline">\Omega </math> the stabilizing effect is more and <math display="inline">\vert {Ra}_2\vert </math> approaching to zero. The result related to linear theory may be compared with the studies of Kiran et al. [[#b0195|[39]]].

==6. Conclusions==

The effect of gravity modulation and vertical throughflow on Bénard–Darcy convection is investigated for linear and weakly nonlinear studies. The porous medium is internally heated. The Ginzburg–Landau model used for weakly nonlinear theory and Venezian approach is used for linear theory. The following conclusions are drawn on previous analysis:
* Effect of <math display="inline">{Pr}_D</math> is to enhance the heat transport for lower values of <math display="inline">\Omega </math> and diminishes for large values of <math display="inline">\Omega </math>.
* Effect of anti-gravity flow ''Pe'' > 0 increases heat transfer and pro-gravity flow ''Pe'' < 0 diminishes heat transfer.
* Upon increasing <math display="inline">\delta </math>, heat transport increases.
* Upon increasing the value of frequency of gravity modulation <math display="inline">\Omega </math>, the amplitude of convection decreases, thus effect of g-jitter becomes negligible at higher values of <math display="inline">\Omega </math>.
* The magnitude of streamlines increases with time and archives maximus size for further values of time.
* Initially isotherms are flat due to conduction state, becomes contour showing the convective regime.
* Low frequency of g-jitter shows its maxim either stabilizing or destabilization effect but, high frequency stabilizes the system.
* Heat source (<math display="inline">R_i>0</math>) and anti-gravity flow (<math display="inline">Pe>0</math>) have destabilizing effect, heat sink (<math display="inline">R_i<0</math>) and pro-gravity flow (<math display="inline">Pe<0</math>) have stabilizing effect on the system.

==Acknowledgements==

The author thank the referees for their valuable suggestions that made the paper into its current form. The author Dr. P. Kiran would like to thank Professor Y. Narasimhulu (Vice Chancellor, RU) for providing research specialities in the department and his encouragement toward research. The author thank Prof. BS. Bhadauria, Ph.D. Supervisor, Institute of Science, Dept of Mathematics, BHU, Varanasi; for motivation toward research.

==References==

<ol style='list-style-type: none;margin-left: 0px;'><li><span id='b0005'></span>
[[#b0005|[1]]] E.S. Nelson; An examination of anticipation of g-jitter on space station and its effects on materials processes; NASA TM (1991), p. 103775</li>
<li><span id='b0010'></span>
[[#b0010|[2]]] M. Wadih, B. Roux; Natural convection in a long vertical cylinder under gravity modulation; J. Fluid Mech., 193 (1988), pp. 391–415</li>
<li><span id='b0015'></span>
[[#b0015|[3]]] H. Wadih, N. Zahibo, B. Roux; Effect of gravity jitter on natural convection in a vertical cylinder; Low-gravity fluid dynamics and transport phenomena (A92-17051 04-29). Washington, DC, American Institute of Aeronautics and Astronautics, Inc. (1990), pp. 309–354</li>
<li><span id='b0020'></span>
[[#b0020|[4]]] G.Z. Gershuni, E.M. Zhukhovitskii, I.S. Iurkov; On convective stability in the presence of periodically varying parameter; J. Appl. Math. Mech., 34 (1970), pp. 442–452</li>
<li><span id='b0025'></span>
[[#b0025|[5]]] P.M. Gresho, R.L. Sani; The effect of gravity modulation on the stability of a heated fluid layer; J. Fluid Mech., 40 (1970), pp. 783–806</li>
<li><span id='b0030'></span>
[[#b0030|[6]]] W.M. Yang; Stability of viscoelastic fluids in a modulated gravitational field; Int. J. Heat Mass Transf., 40 (1997), pp. 1401–1410</li>
<li><span id='b0035'></span>
[[#b0035|[7]]] M.S. Malashetty, V. Padmavathi; Effect of gravity modulation on the onset of convection in a fluid and porous layer; Int. J. Eng. Sci., 35 (1997), pp. 829–839</li>
<li><span id='b0040'></span>
[[#b0040|[8]]] D.A.S. Rees, I. Pop; The effect of G-jitter on vertical free convection boundary-layer flow in porous media; Int. Commun. Heat Mass Transf., 27 (2000), pp. 415–424</li>
<li><span id='b0045'></span>
[[#b0045|[9]]] S. Govender; Weak non-linear analysis of convection in a gravity modulated porous layer; Transp. Porous Media, 60 (2005), pp. 33–42</li>
<li><span id='b0050'></span>
[[#b0050|[10]]] A.V. Kuznetsov; Linear stability analysis of the effect of vertical vibration on bioconvection in a horizontal porous layer of finite depth; J. Porous Media, 9 (2006), pp. 597–608</li>
<li><span id='b0055'></span>
[[#b0055|[11]]] B.S. Bhadauria, P. Kiran; Weak nonlinear double diffusive magneto-convection in a Newtonian liquid under gravity modulation; J. Appl. Fluid Mech., 8 (4) (2014), pp. 735–746</li>
<li><span id='b0060'></span>
[[#b0060|[12]]] B.S. Bhadauria, P.G. Siddheshwar, K. Jogendra, Om.P. Suthar; Nonlinear stability analysis of temperature/gravity modulated Rayleigh–Beńard convection in a porous medium; Transp. Porous Media, 92 (2012), pp. 633–647</li>
<li><span id='b0065'></span>
[[#b0065|[13]]] B.S. Bhadauria, I. Hashim, P.G. Siddheshwar; Study of heat transport in a porous medium under G-jitter and internal heating effects; Transp. Porous Media, 96 (2013), pp. 21–37</li>
<li><span id='b0070'></span>
[[#b0070|[14]]] B. S Bhadauria, P. Kiran; Nonlinear thermal Darcy convection in a nanofluid saturated porous medium under gravity modulation; Adv. Sci. Lett., 20 (2014), pp. 903–910</li>
<li><span id='b0075'></span>
[[#b0075|[15]]] P. Kiran; Nonlinear thermal convection in a viscoelastic nanofluid saturated porous medium under gravity modulation; Ain Shams Eng. J. (2015) [http://dx.doi.org/10.1016/j.asej.2015.06.005 http://dx.doi.org/10.1016/j.asej.2015.06.005]</li>
<li><span id='b0080'></span>
[[#b0080|[16]]] G.M. Homsy, A.E. Sherwood; Convective instabilities in porous media with throughflow; AIChE J., 22 (1976), pp. 168–174</li>
<li><span id='b0085'></span>
[[#b0085|[17]]] M.C. Jones, J.M. Persichetti; Convective instability in packed beds with throughflow; AIChE J., 32 (1986), pp. 1555–1557</li>
<li><span id='b0090'></span>
[[#b0090|[18]]] I.S. Shivakumara, S. Sureshkumar; Convective instabilities in a viscoelastic-fluid-saturated porous medium with throughflow; J. Geophys. Eng., 4 (2007), pp. 104–115</li>
<li><span id='b0095'></span>
[[#b0095|[19]]] I.S. Shivakumara, S. Sureshkumar; Convective instabilities in a viscoelastic-fluid-saturated porous medium with throughflow; J. Geophys. Eng., 4 (2007), pp. 104–115</li>
<li><span id='b0100'></span>
[[#b0100|[20]]] A. Khalili, I.S. Shivakumara; Onset of convection in a porous layer with net throughflow and internal heat generation; Phys. Fluids, 10 (1998), pp. 315–317</li>
<li><span id='b0105'></span>
[[#b0105|[21]]] D.A. Nield; Convective instability in porous media with throughflow; AIChE J., 33 (1987), pp. 1222–1224</li>
<li><span id='b0110'></span>
[[#b0110|[22]]] S.P. Suma, Y.H. Gangadharaiahand, R. Indira; Effect of throughflow and variable gravity field on thermal convection in a porous layer; Int. J. Eng. Sci. Technol., 3 (2011), pp. 7657–7668</li>
<li><span id='b0115'></span>
[[#b0115|[23]]] P. Kiran; Throughflow and g-jitter effects on binary fluid saturated porous medium; Appl. Math Mech., 36 (10) (2015), pp. 1285–1304</li>
<li><span id='b0120'></span>
[[#b0120|[24]]] P. Kiran; Throughflow and gravity modulation effects on heat transport in a porous medium; J. Appl. Fluid Mech., 9 (4) (2015)</li>
<li><span id='b0125'></span>
[[#b0125|[25]]] P. Kiran, B.S. Bhadauria; Nonlinear throughflow effects on thermally modulated porous medium; Ain Shams Eng. J. (2014) [http://dx.doi.org/10.1016/j.asej.2015.03.010 http://dx.doi.org/10.1016/j.asej.2015.03.010]</li>
<li><span id='b0130'></span>
[[#b0130|[26]]] R.K. Vanishree; Effects of throughflow and internal heat generation on a thermo convective instability in an anisotropic porous medium; J. Appl. Fluid Mech., 7 (2014), pp. 581–590</li>
<li><span id='b0135'></span>
[[#b0135|[27]]] G. Hetsroni, M. Gurevich, R. Rozenblit; Natural convection in metal foam strips with internal heat generation; Exp. Therm. Fluid Sci., 32 (2008), pp. 1740–1747</li>
<li><span id='b0140'></span>
[[#b0140|[28]]] B.S. Bhadauria, I. Hashim, P.G. Siddheshwar; Study of heat transport in a porous medium under G-jitter and internal heating effects; Transp. Porous Media, 96 (2013), pp. 21–37</li>
<li><span id='b0145'></span>
[[#b0145|[29]]] B.S. Bhadauria, A. Kumar, J. Kumar, N.C. Sacheti, P. Chandran; Natural convection in a rotating anisotropic porous layer with internal heat; Transp. Porous Media, 90 (2011), pp. 687–705</li>
<li><span id='b0150'></span>
[[#b0150|[30]]] S. Saravanan; Thermal non-equilibrium porous convection with heat generation and density maximum; Transp. Porous Media, 76 (2009), pp. 35–43</li>
<li><span id='b0155'></span>
[[#b0155|[31]]] B.S. Bhadauria, P. Kiran; Heat transport in an anisotropic porous medium saturated with variable viscosity liquid under temperature modulation; Transp. Porous Media, 100 (2013), pp. 279–295</li>
<li><span id='b0160'></span>
[[#b0160|[32]]] W.V.R. Malkus, G. Veronis; Finite amplitude cellular convection; J. Fluid Mech., 4 (1958), pp. 225–260</li>
<li><span id='b0165'></span>
[[#b0165|[33]]] G. Venezian; Effect of modulation on the onset of thermal convection; J. Fluid Mech., 35 (1969), pp. 243–254</li>
<li><span id='b0170'></span>
[[#b0170|[34]]] S.H. Davis; The stability of time periodic flows; Annu. Rev. Fluid Mech., 8 (1976), pp. 57–74</li>
<li><span id='b0175'></span>
[[#b0175|[35]]] A. Srivastava, B.S. Bhadauria, P.G. Siddheshwar, I. Hashim; Heat transport in an anisotropic porous medium saturated with variable viscosity liquid under G-jitter and internal heating effects; Transp. Porous Media, 99 (2013), pp. 359–376</li>
<li><span id='b0180'></span>
[[#b0180|[36]]] R. Kippenhahn, A. Weigert; Stellar Structure and Evolution; Springer (1994)</li>
<li><span id='b0185'></span>
[[#b0185|[37]]] B.S. Bhadauria, P. Kiran; Weak nonlinear analysis of magnetoconvection under magnetic field modulation; Phys. Scripta, 89 (9) (2014), p. 095209</li>
<li><span id='b0190'></span>
[[#b0190|[38]]] B.S. Bhadauria, P. Kiran; Effect of rotational speed modulation on heat transport in a fluid layer with temperature dependent viscosity and internal heat source; Ain Shams Eng. J., 5 (4) (2014), pp. 1287–1297</li>
<li><span id='b0195'></span>
[[#b0195|[39]]] P. Kiran, B.S. Bhadauria, V. Kumar; Thermal convection in a nanofluid saturated porous medium with internal heating and gravity modulation; J. Nanofluid (2016) accepted</li>
</ol>