Nonlinear stability analysis of a fluid saturated anisotropic Darcy–Brinkman medium with internal heat source

Applied Mathematics and Computation 358 (2019) 216–231

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Nonlinear stability analysis of a fluid saturated anisotropic Darcy–Brinkman medium with internal heat source Reena Nandal∗, Amit Mahajan Department of Applied Sciences, National Institutes of Technology, Narela, Delhi 110040, India

a r t i c l e

i n f o

Keywords: Anisotropy Internal heat parameter Darcy–Brinkman model Energy method

a b s t r a c t In the present study, the effect of mechanical anisotropy on the onset of convection is observed when the fluid layer is influenced by internal heat generation and heat extraction through the boundary. The extended Darcy model is employed to establish the momentum equation. Particularly, the emphasis is laid upon that how these factors alter the criteria of onset of convection. The variation in certain parameters affects the stability of the system. Linear instability and nonlinear stability analysis are performed to predict the stability of the system, however; it is found that the linear and nonlinear analyses are not in agreement with the substantial change in internal heat generation. All numerical and graphical results of the subsequent analysis are obtained by using the Chebyshev-pseudospectral method. © 2019 Elsevier Inc. All rights reserved.

1. Introduction Thermal convection in a plane layer in the presence of an internal heat source is a well-defined phenomenon. The presence of heat source leads to the phenomenon of penetrative convection where part of the layer has a tendency to move whereas the remainder of the layer will remain motionless until a certain point when movement in the rest of the layer penetrates into the stable layer and a resultant motion then ensues. Penetrative convection has immense application in geophysics, atmospheric physics, building design, biochemical decay, it occurs in Sun and responsible for volcanic plumes into the Earth’s atmosphere [Straughan [49] and references therein]. In a thermodynamical system, internal heat generation studies are performed most extensively in the context of two mechanisms; one is a nonlinear density temperature relationship which accommodates the maximum density character of fluid and second is a heat source or sink. Natural convection due to internal heat generation in natural and engineering phenomenon can happen for many reasons. For example, it is observed for a laboratory model, an electrolyte can cause internal heat generation, internal energy generation due to absorption of sunlight, is responsible for cellular convection in an atmosphere of Venus and radioactive decay of isotopes causes the earth’s mantle convection [Kuznetsov and Nield [26] and references therein], biochemical decay causes internal heat generation [36], electronic thermal conductivity is responsible for deep mantle convection of exosolar planets [52], rise of volcanic plumes due to penetrative convection into the earth’s atmosphere [24]. Roberts [40] has made a significant contribution to the study of convection in a horizontal fluid layer which is cooled from above, thermally insulated from below and heated uniformly by an internal heat source. Very sparse literature is available for the study of penetrative convection and a huge qualitative analysis is provided by researchers with varying mathematical models. Matthews [31], Straughan [48,47],

∗

Corresponding author. E-mail address: [email protected] (R. Nandal).

https://doi.org/10.1016/j.amc.2019.03.023 0 096-30 03/© 2019 Elsevier Inc. All rights reserved.

R. Nandal and A. Mahajan / Applied Mathematics and Computation 358 (2019) 216–231

217

Ames and Cobb [2], Srimani and Sudhakar [44], Herron [23], Carr and Straughan [11], Berlengiero [3], Goluskin and Spiegel [18], Hattori et al. [22], Goluskin et al. [19], Mahajan and Sharma [27] are the few who provides the well documentation of the topic and observe the penetrative convection at distinct aspects. Recently, Nandal and Mahajan [33] examined the effect of penetrative convection in a Darcy–Brinkman porous medium when the fluid layer is influenced by four internal heat generating functions. Convection in porous media is generally dealt with isotropy but there are many situations where anisotropy exists in the medium and these situations are encountered in numerous systems in industry and nature. If in the porous media the fluid permeates non-uniformly due to the direction of fluid flow or pressure gradient then the porous medium is considered as anisotropic. The variation in permeability reflects that there is the difference in path length by which fluid particle moves and the forces it experiences during movement. Many porous media encountered in nature and industry are not isotropic e.g. rock, sand, wood, sponge, fibrous insulating materials, filter bed, river bed, oil sand, underground rock formation etc. The measurement of macroscopic transport property of a fluid in porous media such as permeability, diffusivity, and thermal conductivity are direction dependent, so the anisotropy of the porous medium highly influences the property of fluid and flow pattern in a porous medium. Anisotropy is generally a consequence of preferential orientation and asymmetric geometry of the grain which constitute the porous media and it is very difficult to characterize it qualitatively [39]. Castinel and Combarnous [12] first studied in this area and obtained the criterion for the onset of convection in a fluid layer saturated in anisotropic porous medium with impermeable lower and upper boundaries. The anisotropy in the porous medium can be seen by two facts that are by variation in directional permeability or by variation in directional thermal conductivity or by both. Castinel and Combarnous [12] and Epherre [15] commented that if the vertical permeability and vertical thermal conductivity keep fixed and increase the horizontal permeability, then variation in horizontal permeability promotes the horizontal movement which increases the cell width and critical Rayleigh number get reduces. On the other hand increase in horizontal thermal conductivity decays the buoyancy forces which results in the increase in cell width and critical Rayleigh number. The pioneer work of Tyvand and Storesletten [51] is a very beautiful paper which provides a deep knowledge of the topic anisotropy. Tyvand and Storesletten observed that in the isotropic situation, the rolls of preferred motion are with a square cross section but for the anisotropic medium the preferred motion in tilted cells. Further, Storesletten [45] consider the anisotropic thermal conductivity for the observation of free convection in a porous medium. Straughan [46] used the same configuration of Tyvand and Storesletten [51] but for a non-Boussinesq (nonlinear temperature density function) model to observe the penetrative convection. The onset of Rayleigh Bénard convection is investigated by Qin and Kaloni [37] in a porous medium with anisotropic permeability. Various researchers [Mahidjiba et al. [29], Malashetty and Kollur [30], Guo [20], Gaikwad and Dhanraj [16], Safi and Benissaad [41], Karmakar and Sekhar [25]] perform their research in this area and provide a better platform to the learner to understand the phenomenon of anisotropy. Penetrative convection in anisotropic porous media is studied considerably because of its occurrence in various geophysical and industrial processes. A good amount of work can be obtained on this topic in literature. The combined effect of internal heat generation and inclined temperature gradient on the thermal instability is observed by Parthiban and Patil [35] in an anisotropic porous media. Capone et al. [10] studied the penetrative convection in a porous medium for the internal heat source and sink for constant anisotropic thermal diffusion and non-homogeneous permeability. The effects of variable permeability and thermal diffusivity on the onset of anisotropic penetrative convection are observed by Capone et al. [9], for fluids obeying quadratic density law. A study is taken by Bhadauria [4] on double diffusive convection in an anisotropic porous media when convection is induced by internal heat source and fluid heated and salted from below with generalized Darcy model employed on momentum equation. Harfash [21] did the stability analysis for penetrative convection in an anisotropic porous medium for the fluid obeyed quadratic density law with anisotropic permeability and thermal diffusivity effects. The effect of internal heat generation on the onset of convection in a nanofluid saturated anisotropic porous media is investigated by Shivakumara and Dhananjaya [42]. Degan et al. [14] observed the effects of hydrodynamic anisotropy on the mixed convection in a vertical Darcy–Brinkman porous channel heated on its plate with thermal radiations for a fully developed flow regime. The work of Capone et al. [8], Shivakumara et al. [43], Altawallbeh et al. [1], Gangadharaiah et al. [17], Mahajan and Nandal [28] contributes a lot to understand the topic penetrative convection in anisotropic porous media. In the present paper, we observed the effects of the internal heat source and heat extraction through the boundary on the onset of Darcy–Brinkman [6] convection in an anisotropic porous medium. The value of critical Rayleigh number for linear and nonlinear analysis is found very close when penetration is effective near the middle part of the fluid layer. Effects of anisotropy of the porous medium, Darcy–Brinkman number, internal heat source parameter, on the critical value of Rayleigh number by linear and nonlinear analysis are obtained numerically and graphically by Chebyshev tau-QZ method. 2. Mathematical formulation of problem The physical configuration consists of a thin layer of an incompressible fluid of thickness d in an anisotropic porous medium with an internal heat source Q(z). The fluid is confined between z ∈ (0, d) with ∂∂Tz = γ at z = 0 and T = TU at z = d, here γ represents the rate at which heat is extracted through the lower boundary. Gravity g acts vertically downwards in the negative z-direction and we considered that the density has a linear dependence on temperature as:

ρ f (T ) = ρ0 [1 − α (T − T0 )]

(1)

218

R. Nandal and A. Mahajan / Applied Mathematics and Computation 358 (2019) 216–231

Fig. 1. Internal heat source profile with the depth (z) of fluid layer.

where ρ 0 and T0 is the reference density and temperature and α is the thermal expansion coefficient. Hence the governing equations for the mechanically anisotropic porous medium are as follows: Momentum Balance:

ρf ∂q ˆ − M i j μq + μ ∇ 2 q , = −∇ P + ρ f (T )gk ε ∂t

(2)

Mass Balance:

∇ · q = 0,

(3)

Temperature Equation:

A

∂T + q · ∇ T = κ∇ 2 T + Q. ∂t

(4)

where ρ f , ε , q, P, μ, μ , t, A and κ are the density of the fluid, porosity, filter velocity, pressure, dynamic viscosity, effective viscosity, time, ratio of heat capacity and thermal diffusivity, respectively and Mij is inverse tensor to permeability     K1 0 0 ζ 0 0 1 tensor Kij , that is MK = I with Ki j = 0 K1 0 and Mi j = K 0 ζ 0 where ζ = KK1 [Tyvand and Storesletten [51], 0 0 K 0 0 1 Straughan [50]]; which represents the anisotropy parameter and K1 , K are horizontal permeability and vertical permeability of the medium. The term Q represents the internal heat source and its presence induces the penetrative convection in the porous layer. In a system, heat can be generated by many ways. Here we considered four types of internal heat generating functions that are constant, linearly increasing, decreasing and non-uniform. The constant heat generating functions is considered by various authors in the early studies of penetrative convection for a laboratory model. But in general, there are many processes such as biochemical decay, microwave heating, earth’s mental convection, nuclear reactor processes, chemical reactions and many more geophysical and engineering processes, the internal heat generation is not in a constant way. So, we considered the other type of internal heat generating functions that may occur in the above discussed physical problems. Other examples of penetrative convection induced by internal heating can be found in Roberts [40], Matthews [32] and Straughan [47]. So the following internal heat generating functions are considered to observe the effect of internal heat generation on the onset of convection in porous media. Case A: Q = Q0 .

1



z . 2 d   3z 2 z Case C: Q = Q0 2 + − 3 . d 2d 2   2π z 4π z Case D: Q = Q0 1 + sin + sin . d d Case B: Q = Q0

+

The dimensionless profile of these heat generating functions with the depth of fluid layer is displayed in Fig. 1. Now, the solution to the governing equations in quiescent state is given by

qb = 0, Pb = P (z ), ρb = ρ (z ), Tb = T (z ),

(5)

R. Nandal and A. Mahajan / Applied Mathematics and Computation 358 (2019) 216–231

219

where subscript ‘b’ represents the basic state and basic temperature function Tb is given as:

Tb =

⎧ Q0 2 ⎪ − (z − d2 ) + γ z + TU − γ d, ⎪ ⎪ 2 κ ⎪  ⎪ ⎪ Q0 z 3 5 2 z2 ⎪ ⎪ ⎪ ⎨− κ 6d + 4 − 12 d + γ z + TU − γ d,   4

case A case B (6)

3

Q0 z 5 z ⎪ − − + z2 − d2 + γ z + TU − γ d, ⎪ 2 ⎪ κ 2 d 8 8 d ⎪ ⎪  2  ⎪ 2 ⎪ Q d d2 d2 3d 2 z 2 π z 4π z 3d ⎪ 0 ⎪ − sin − sin + z− − + γ z + TU − γ d. ⎩− κ 2 4π 2 d d 4π 2 4π 16π 2

case C case D

For the investigation of the effect of penetrative convection, we have to obtain a certain range of a parameter i.e. internal heat parameter ξ [defined later]. The range can be obtained from the extremum of density. The obtained range of ξ for all type of internal heat generating function is (0, 1) and outside this range, convection is purely due to buoyancy forces. Fig. 2 represents the profile of basic dimensionless temperature (Tb ) with depth of fluid layer (z) for heat input through fluid layer and heat extracted through the boundary. Now to analyze the stability of the system, the following perturbations are introduced in the steady state solution:

q = qb + q , P = Pb + P  , T = Tb + θ , ρ = ρb + ρ  ,

(7)

where q , P , θ , ρ  are the perturbed quantity in velocity, pressure, temperature, and density. Now, the following nonlinear system of perturbation equations (after dropping prime) is obtained:

ρf ∂q ˆ − M i j μq + μ ∇ 2 q , = −∇ P + ρ f gαθ k ε ∂t

(8)

∇ · q = 0,

(9)

A

dT ∂θ + q · ∇θ = − b w + κ∇ 2 θ . ∂t dz

(10)

To reduce the complexity of the above prototype system it is convenient to write the above system of equation in dimensionless form. Accordingly, we introduced the following scaling parameter:

z = d z∗ , t =

d2 A

κ

t ∗, q =

κ d

q∗ , P =

μκ K

P∗ ,

θ=

⎫

Q0 d 2 ∗ ⎪ θ ,⎪ κ Ra ⎬

ρ f α g|Q0 |d3 K ˜ εν d2 A K μ K ⎪ ⎪ ⎭ Va = , Da = 2 , Ra = , Da = 2 κK d μκ μd 2

(11)

where Va is Vadasz number, Da is Darcy number, Ra is a thermal Darcy-Rayleigh number, D˜ a is Darcy–Brinkman number. Let mij = diag(ζ , ζ , 1) and hence the system after non-dimensionalization leads to the following nondimensional perturbed equations (after dropping ∗ )

1 ∂q ˆ, = −∇ P − mi j q + D˜ a∇ 2 q + (Ra )1/2 θ k Va ∂t

(12)

∇ · q = 0,

(13)

∂θ + q · ∇ θ = ∇ 2 θ + m(z )(Ra )1/2 w. ∂t

(14)

where the function m(z) represents non-dimensional temperature gradient functions for case A-D, its value is:

⎧ ⎪ ⎪HH1 z− H2 ξ, ⎪ 1 ⎪ ⎪ ⎨ 2 z2 + z − H2 ξ , 1  3 m (z ) = H1 z3 − z2 + 2z − H2 ξ , ⎪ ⎪ 2 ⎪  2 ⎪ ⎪ ⎩H1 z − 1 cos 2π z − 1 cos 4π z + 3 − H2 ξ . 2π 4π 4π

case A case B case C

(15)

case D

where H1 = sign(Q0 ), H2 = sign(γ ), represents the heat generation within the fluid layer and heat extraction through the κ|γ | lower boundary respectively. ξ = |Q |d is internal heat parameter which represents the strength of heat due to heat gener0

ation and extraction. Eqs. (12)–(14) are to be solved for the following impermeable and isothermal boundary conditions [Straughan [47]]

q = 0 at z = 0, 1.

(16)

θ = 0 at z = 1 and Dθ = 0 at z = 0.

(17)

220

R. Nandal and A. Mahajan / Applied Mathematics and Computation 358 (2019) 216–231

Fig. 2. Basic temperature profiles (Tb ) with depth of fluid layer (z) when heat input through the fluid layer (H1 = 1) and heat extracted through the lower boundary (H2 = 1).

R. Nandal and A. Mahajan / Applied Mathematics and Computation 358 (2019) 216–231

221

3. Linear analysis The main impetus behind the linear analysis is that the linear theory provides boundaries for instability and the potential growth rate of perturbations [Chandrasekhar [13]]. For linear analysis, we return to the perturbed dimensionless Eqs. (12)– normal mode analysis assuming normal modes of (14) by discarding nonlinear terms. We perform the standard stationary  the form (q, θ ) = [q(z), (z)]exp [σ t + i(ax x + ay y)], where a = (a2x + a2y ) is the wave number of disturbance and σ is a d function of wave number and other parameters of the system. Letting D ≡ dz and to eliminate the pressure term, operate double curl on the equation (11). By considering the third component of the resulting equation, we have the following linearized system of equations:

−

σ  Va





2

D2 − a2 W = −D˜ a D2 − a2 W +



 ζ D2 − a2 W + a2 (Ra )1/2 ,

  σ = m(z )(Ra )1/2W + D2 − a2 .

(18)

(19)

For numerical analysis, we use the Chebyshev pseudospectral-QZ method. The domain of method is [ − 1, 1], so we reset the present domain from [0, 1] with coordinate transformation z to 2z − 1, hence the resulting equations are:

−

σ  Va





2





4D2 − a2 W = −D˜ a 4D2 − a2 W + 4ζ D2 − a2 W + a2 (Ra )1/2 ,

  σ = n(z )(Ra )1/2W + 4D2 − a2 .

(20)

(21)

where n(z) for case A–D is given as

⎧ z + 1 ⎪ H1 − H2 ξ , ⎪ ⎪ ⎪  z 2+ 1  z + 3  ⎪ H ⎪ ⎪ − H2 ξ , ⎨ 1 2 2 2     n (z ) = 3 z+1 2 1 z+1 3 ⎪ H z + 1 − + − H2 ξ , ( ) ⎪ 1 ⎪ 2 2 2 2 ⎪ ⎪     ⎪ ⎪ ⎩H1 z + 1 + 1 cos π z − 1 cos 2π z + 3 − H2 ξ . 2 2π 4π 4π

case A case B (22) case C case D

and the corresponding impermeable, isothermal boundary leads to the following conditions:

W = 0, DW = 0 at z = −1, 1.

(23)

= 0 at z = 1 and D = 0 at z = −1.

(24)

For the solution of Eqs. (20) and (21) with boundary conditions (23) and (24), we expand the unknown function W and

in Chebyshev polynomials Tn (z) = cos (ncos −1 z) on [ − 1, 1] as [Canuto et al. [7]]: W (z ) =

L 

bn Tn (z )

(25)

n=0

where

bn =



2 cn π

1

−1

W (z )Tn (z ) dz √ 1 − z2

(26)

here c0 = 2 and ci = 1(i ≥ 1). Now we replace the integration in (26) with numerical discrete integral on the collocation points,

bn =

L 2  1     W z j Tn z j cn L cj

(27)

j=0

where zj represents collocation points and z j = − cos( jπ /L ) for j = 1, 2, 3... . Now Eq. (25) can be expressed as:

W (z ) =

L 

 

m j (z )W z j

j=0

where the interpolating polynomial mj (z) is given as

(28)

222

R. Nandal and A. Mahajan / Applied Mathematics and Computation 358 (2019) 216–231

Fig. 3. Neutral Stability curve at distinct value of internal heat parameter (ξ ) and anisotropy parameter (ζ ) for case A.

Fig. 4. Neutral Stability curve at distinct value of internal heat parameter (ξ ) and anisotropy parameter (ζ ) for case B.

m j (z ) =

L 2  1   Tn z j Tn (z ) Lc j cn

(29)

n=0

m j  ( zi ) =

L 2  1    T z j Tk (zi ) c jL ck k

(30)

k=0

the Eq. (30) represents the differentiation matrix D and for D2 ,D3 ,D4 a relationship is of the type Dn + 1 = DDn . Applying similar calculation for and after appropriate discretization Eqs. (20) and (21) with boundary conditions (23) and (24), can be written in matrix form as

σ BX = A1 X,

(31)

where B and A1 are complex matrices, given as:



B=



1 − Va 4 D2 − a2 I O





O , A1 = I





2



−D˜ a 4D2 − a2 I + 4ζ D2 − a2 I n(Z )(Ra )1/2



 1/2 2 a (2Ra ) 2 I . 4D − a I

R. Nandal and A. Mahajan / Applied Mathematics and Computation 358 (2019) 216–231

223

Fig. 5. Neutral Stability curve at distinct value of internal heat parameter (ξ ) and anisotropy parameter (ζ ) for case C.

Fig. 6. Neutral Stability curve at distinct value of internal heat parameter (ξ ) and anisotropy parameter (ζ ) for case D.

where O and I are zero and identity matrices. Numerical results for the linear theory are presented in Section 5. 4. Nonlinear stability analysis For the nonlinear stability analysis of the system, the energy method is one of the oldest methods. Often it is related to the various types of energy measurements that are why it named as energy method. The linear analysis often provides little information on the behavior of the nonlinear system as any potential growth in nonlinear terms is not considered. In order to analyze the behavior of the nonlinear system, one considers finite perturbations and retains the quadratic terms. The system may be analyzed by generalized nonlinear energy technique to obtain stability threshold which may be compared with instability threshold obtained by linear theory. This method can trace to the work of Reynolds [38] and Orr [34] and one can surely achieve the condition under which energy will decrease. In energy method, energy can be constructed by using coupling parameter [Joseph (1966)]. The energy method is found to be successful in a variety of contexts such as electrodynamics, magneto-hydrodynamics, convection with chemical reactions, ferrohydrodynamics, half-space problems, surface tension drove convection problems, geophysical problems, thermohaline circulation in oceans, multi-component convection, and mathematical biology. Energy method seeks to find the configurations of minimum potential energy and then attempt to

224

R. Nandal and A. Mahajan / Applied Mathematics and Computation 358 (2019) 216–231

Fig. 7. Variation of critical Rayleigh Number (Ra) with internal heat parameter (ξ ), for case A of constant heat source (H1 = 1) and heat extracted through the lower boundary (H2 = 1) by linear and nonlinear analysis.

establish it as the stable configuration. Nonlinear stability energy method has been employed highly successfully in findings sharp bounds in hydrodynamical problems. Also, it is useful to delimit the parameter region of possible subcritical instability. Virtue of this method is that it gives insight into the physical situation and it can accommodate for finite disturbances. Due to the numerous and diverse application, the energy method is undoubtedly important. Now, in order to analyze the nonlinear stability, Eq. (12) is multiplied by q and Eq. (14) with θ and integrated over V (after using Eq. (13), divergence theorem and boundary conditions) to obtain





d (Ra )1/2 I E (t ) = −D 1 − . dt D

(32)

λ q 2 , and I = λ wθ  + m(z)wθ , is energy production term and D = λm q 2 + ∇θ 2 + where E (t ) = 12 θ 2 + 2Va ij 2 D˜ aλ ∇ q is the energy dissipating terms. where • denotes the L2 (V) norm, •  denotes the integration and V is the typical periodicity cell. Next, define RaE by

1 I = max , RaE  D

(33)

where  is the space of admissible solutions, such that  = {q, θ |q, θ ∈ L2 [0,1]}, in addition, q and θ satisfies the boundary condition (16) and (17).

R. Nandal and A. Mahajan / Applied Mathematics and Computation 358 (2019) 216–231

225

Fig. 8. Variation of critical Rayleigh Number (Ra) with internal heat parameter (ξ ), for case B of increasing heat source (H1 = 1) and heat extracted through the lower boundary (H2 = 1) by linear and nonlinear analysis.

Now, following the Eqs. (32) and (33), we may write



d (Ra )1/2 E (t ) ≤ −D 1 − dt RaE



,

(34)

If (Ra )1/2 < RaE , using Poincaré’s inequality in (34) leads to the exponential decay of energy with time. Now, RaE is to be determined from (33) to find the condition of nonlinear stability. So, to find the maximum of (33) we use the calculus of variation and the maximizing solution satisfies the associated Euler-Lagrange equation, given as:

where

δ D − RaE δ I = 0,

(35)

    δ D = 2λmi j q · h − 2 ∇ 2 θ η − 2D˜ aλ ∇ 2 q · h .

(36)

δ I = λ wη + θ h3  + m(z )(wη + θ h3 ).

(37)

for all admissible h and η. Since h and η are arbitrary, the Eq. (35) is equivalent to

ˆ − RaE m(z )θ k ˆ = 2∇ p, 2λmi j q − 2D˜ aλ∇ 2 q − RaE λθ k

(38)

226

R. Nandal and A. Mahajan / Applied Mathematics and Computation 358 (2019) 216–231 Table 1 Comparison of the results with the existing studies. Internal Heat Parameter

Present Study

Existing Result

ξ = 0; D˜ a = 0 ξ = 0; D˜ a = 0.2

RaL = 471.385 (for case A) RaL = 3994.759 (for case A)

RaL = 471.38466 [Borujerdi et al. [5]] RaL = 3994.761 [Nandal and Mahajan [ 33 ]]

Table 2 The variation of Rayleigh Number (Ra) with internal heat parameter (Q) at fixed values of Darcy–Brinkman number D˜ a = 0 and different values of anisotropy parameter (ζ ) for Q - source and γ - sink by linear/nonlinear analysis for case A, B, C, and D. D˜ a = 0

ζ

A

ξ

B

C

D

acL

RaL

acE

λ

RaE

acL

RaL

acE

λ

RaE

acL

RaL

acE

λ

RaE

acL

RaL

acE

λ

RaE

0.125 0.1 0.2 0.3 0.4 0.5 0.6

1.38 1.41 1.56 2.17 2.80 3.50

28.06 40.04 65.60 118.07 209.90 410.16

1.39 1.43 1.52 1.74 2.07 2.15

0.36 0.28 0.22 0.18 0.17 0.22

26.70 36.59 54.92 91.28 161.89 276.37

1.41 1.56 2.24 2.92 3.64 4.71

40.52 67.22 121.71 210.26 383.50 786.28

1.43 1.53 1.79 2.18 2.34 2.29

0.27 0.21 0.18 0.18 0.20 0.29

37.19 56.33 94.15 165.15 284.69 431.23

1.36 1.37 1.38 1.44 1.71 2.35

18.83 23.79 32.09 48.20 84.52 159.18

1.37 1.38 1.40 1.44 1.56 1.79

0.50 0.41 0.33 0.25 0.19 0.16

18.37 22.90 30.10 42.75 67.63 120.04

1.36 1.36 1.37 1.38 1.42 1.59

14.56 17.38 21.50 28.05 39.69 63.23

1.37 1.37 1.38 1.40 1.43 1.52

0.64 0.55 0.46 0.37 0.30 0.23

14.26 16.88 20.61 26.23 35.42 51.78

0.25

0.1 0.2 0.3 0.4 0.5 0.6

1.67 1.72 1.93 2.62 3.33 4.16

37.62 53.09 84.93 148.39 261.87 511.38

1.68 1.73 1.85 2.16 2.60 2.77

0.37 0.29 0.22 0.18 0.17 0.20

35.97 48.99 72.93 120.20 214.16 383.24

1.72 1.93 2.70 3.47 4.33 5.60

53.84 87.20 152.97 262.43 478.33 980.61

1.73 1.86 2.18 2.80 3.09 2.97

0.28 0.22 0.19 0.17 0.18 0.27

49.83 74.79 123.59 215.92 384.90 606.31

1.64 1.65 1.68 1.77 2.11 2.81

25.44 31.99 42.81 63.34 107.85 199.02

1.65 1.66 1.69 1.76 1.91 2.21

0.51 0.42 0.33 0.25 0.19 0.16

24.88 30.91 40.45 57.11 89.76 159.99

1.64 1.64 1.65 1.68 1.75 1.97

19.71 23.44 28.87 37.39 52.24 81.33

1.64 1.65 1.66 1.69 1.74 1.84

0.65 0.56 0.47 0.38 0.30 0.24

19.34 22.84 27.80 35.27 47.39 68.95

0.5

0.1 0.2 0.3 0.4 0.5 0.6

2.03 2.12 2.41 3.16 3.96 4.94

53.24 74.20 115.56 195.62 342.84 669.19

2.05 2.12 2.29 2.67 3.33 3.60

0.37 0.29 0.23 0.19 0.16 0.17

51.09 69.04 101.51 164.58 290.99 546.95

2.11 2.41 3.26 4.13 5.14 6.65

75.41 118.81 201.56 343.64 626.10 1283.45

2.12 2.30 2.75 3.43 4.20 4.09

0.29 0.23 0.19 0.18 0.16 0.25

70.28 104.03 168.69 290.93 530.57 879.74

1.99 2.01 2.06 2.19 2.61 3.36

36.32 45.43 60.29 87.71 144.57 261.12

1.99 2.02 2.06 2.15 2.36 2.77

0.52 0.42 0.34 0.26 0.20 0.16

35.58 44.05 57.32 80.24 124.58 220.14

1.98 1.99 2.01 2.06 2.16 2.45

28.18 33.41 40.95 52.61 72.50 110.02

1.98 2.00 2.02 2.06 2.13 2.28

0.66 0.56 0.47 0.39 0.31 0.24

27.70 32.64 39.59 49.98 66.70 96.17

0.75

0.1 0.2 0.3 0.4 0.5 0.6

2.30 2.41 2.75 3.53 4.39 5.47

66.82 92.36 141.47 235.30 410.93 802.00

2.30 2.40 2.60 3.06 3.78 4.37

0.38 0.30 0.24 0.19 0.16 0.16

64.23 86.32 125.74 201.30 352.94 679.49

2.40 2.76 3.64 4.57 5.69 7.36

93.99 145.52 242.36 411.91 750.41 1538.25

2.39 2.63 3.10 3.88 4.93 4.97

0.30 0.23 0.20 0.18 0.15 0.24

87.89 128.81 206.01 352.04 644.36 110 0.0 0

2.23 2.26 2.33 2.49 2.96 3.73

45.84 57.15 75.43 108.57 175.57 313.35

2.24 2.27 2.33 2.44 2.69 3.16

0.52 0.43 0.34 0.27 0.20 0.16

44.94 55.50 71.94 100.08 153.87 269.27

2.22 2.24 2.27 2.33 2.46 2.78

35.61 42.13 51.47 65.81 89.89 134.32

2.22 2.24 2.27 2.32 2.41 2.58

0.66 0.57 0.48 0.39 0.32 0.25

35.03 41.21 49.87 62.74 83.30 119.24

2∇ 2 θ + λRaE w + m(z )RaE w = 0.

(39)

where p is a Lagrange’s multiplier introduced, since q is solenoidal. On taking curlcurl of Eq. (38) to remove the Lagrange’s multiplier, we have the following resulting equation

  ∂ 2w ∇12 w + ζ 2 + 2D˜ aλ∇ 4 w + λRaE ∇12 θ + RaE m(z )∇12 θ = 0. ∂z

−2λ

(40)

Now we assume a plane tiling form (w, θ ) = [W(z), (z)]h(x, y), where h(x, y) is the normalized function which tiles the plane and satisfying ∇12 h + a2 h = 0, a being the wave number (Chandrasekhar [13], pp. 106–114). The wave number is found a posteriori to be non-zero, so from Eqs. (39) and (40), we see that W, satisfy



−2λ



  2 ζ D2 − a2 W + 2D˜ aλ D2 − a2 W − λRaE a2 − a2 m(z )RaE = 0,

(41)



2 D2 − a2 + λRaE W + m(z )RaE W = 0.

(42)

and the corresponding boundary conditions are

W = 0, DW = 0 at z = 0, 1.

(43)

= 0 at z = 1 and D = 0 at z = 0.

(44)

Now, following the method explained for linear analysis the domain is reset from [0, 1] to [ − 1, 1] to match the domain of Chebyshev pseudospectral-QZ method. The system of Eqs. (41) and (44) forms a sixth-order eigenvalue problem which is utilized to locate the nonlinear critical Darcy-Rayleigh number RaE , which is given by RaE = max min RaE (λ, a ). The derived results of energy approach are presented in Section 5.

λ

a

R. Nandal and A. Mahajan / Applied Mathematics and Computation 358 (2019) 216–231

227

Fig. 9. Variation of critical Rayleigh Number (Ra) with internal heat parameter (ξ ), for case C of decreasing heat source (H1 = 1) and heat extracted through the lower boundary (H2 = 1) by linear and nonlinear analysis.

5. Result and discussion The linear and nonlinear results here reported are in a full range data that is extracted from the eigenvalue problem (20)–(24) and (41)–(44). Chebyshev pseudospectral method scheme is adopted to solve the eigenvalue problem. Here we are showing the results for the different types of internal heat generating functions including the constant heat source, increasing, decreasing and non-uniform heat sources with boundary effects. On the theoretical ground as we reduce our system to the case of pure internal heat generation situation when heat is not extracted through the lower boundary and for an isotropic porous media, the results are found in agreement with the existing studies (Table 1). Numerical results relative to the effects of internal heat parameter, Darcy–Brinkman number and anisotropy parameter on the critical value of Rayleigh number are presented in Tables 2–4. Our main impetus during this analysis is that how variation in distinct parameter affects the stability of the system and is there any agreement between the two theories that we used for the stability analysis. We first consider the case A of the constant heat source. It is shown from Table 2 when D˜ a = 0 and ζ = 0.125 (ratio of vertical permeability to the horizontal permeability), with the increase of internal heat parameter value of critical Rayleigh number increases. The reason for this behavior can be understood from Fig. 2. For ξ = 0.1, there is a high-temperature difference between the lower and upper boundary which enhances the rate of movement of fluid particles due to thermal transfer from bottom to top and results in fast convection. But as the strength of ξ increases, it lowers the temperature difference between the boundaries and so the fluid layer becomes more stabilized. Because the

228

R. Nandal and A. Mahajan / Applied Mathematics and Computation 358 (2019) 216–231

Fig. 10. Variation of critical Rayleigh Number (Ra) with internal heat parameter (ξ ), for case D of non-uniform heat source (H1 = 1) and heat extracted through the lower boundary (H2 = 1) by linear and nonlinear analysis.

effect of heat generation within the fluid layer is counterbalanced by the heat extraction through the lower boundary and hence higher number of Rayleigh is required to achieve the onset of convection. In our study, we are considering that 0 < ζ < 1, means horizontal permeability is higher than the vertical permeability. In this situation, a fluid flow easily along the principal axis which is of most physical interest because of it creates a larger pressure drop near the surface of the wellbore. So if the horizontal permeability is high it is very hard to extract the oil or water from the earth. But in account of convective heat transfer as the fluid flow easily along the principal axis the rate of heat transfer will be very fast but as the value of horizontal permeability reduces or anisotropy parameter increases (horizontal permeability approaches vertical permeability), the value of critical Rayleigh number also increases and the system becomes more stable as the results in Tables 3 and 4 reveals. Further, we increase the value of Darcy–Brinkman number, critical Rayleigh number tends to high values which result in the highly stable system. The fluid dynamical system behaves in the same manner for other heat generating functions B, C and D. Wave number increases with the increase of internal heat parameter that is cell size decreases with increase in internal heat parameter. It is further observed that coupling parameter λ does not depend on anisotropy parameter and Darcy–Brinkman number but it depends on the change of strength of internal heat parameter. Figs. 3–6, shows a selection of neutral stability curves for all the four heat generating functions at varying values of internal heat parameter. These figure composed of two frames: one for ζ = 0.25 and other for ζ = 0.75 for Darcy regime (D˜ a = 0). Neutral stability provides the stability criterion in terms of critical Rayleigh number, above the linear plot, the sys-

R. Nandal and A. Mahajan / Applied Mathematics and Computation 358 (2019) 216–231

229

Table 3 The variation of Rayleigh Number (Ra) with internal heat parameter (Q) at fixed values of Darcy–Brinkman number D˜ a = 0.2 and different values of anisotropy parameter (ζ ) for Q - source and γ - sink by linear/nonlinear analysis for case A, B, C, and D. D˜ a = 0.2 A

ζ

ξ

B

C

D

acL

RaL

acE

λ

RaE

acL

RaL

acE

λ

RaE

acL

RaL

acE

λ

RaE

acL

RaL

acE

λ

RaE

0.125 0.1 0.2 0.3 0.4 0.5 0.6

2.59 2.64 2.76 3.17 4.10 5.07

744.33 1014.74 1571.62 3106.45 7982.17 23,705.19

2.60 2.65 2.76 3.01 3.56 3.91

0.39 0.30 0.21 0.15 0.11 0.14

729.39 978.63 1454.86 2562.23 5681.23 13,798.75

2.64 2.77 3.24 4.24 5.24 6.80

1052.53 1662.07 3406.76 8723.68 24,073.35 81,548.36

2.65 2.77 3.08 3.67 4.01 3.87

0.29 0.21 0.14 0.11 0.14 0.23

1013.61 1528.77 2746.60 6093.93 14,397.65 24,970.12

2.55 2.57 2.60 2.66 2.82 3.40

513.32 630.73 816.62 1152.59 1913.59 4350.05

2.57 2.58 2.60 2.67 2.80 3.13

0.40 0.40 0.36 0.26 0.18 0.12

496.44 620.18 798.95 1106.41 1743.33 3444.40

2.54 2.55 2.57 2.60 2.64 2.75

399.02 466.70 561.75 704.67 942.37 1406.45

2.55 2.56 2.58 2.60 2.65 2.74

0.70 0.60 0.50 0.41 0.32 0.23

396.27 462.34 554.24 690.28 909.83 1312.20

0.25 0.1 0.2 0.3 0.4 0.5 0.6

2.60 2.65 2.77 3.18 4.11 5.08

751.73 1024.68 1586.39 3130.78 8022.33 23,785.51

2.61 2.66 2.77 3.03 3.57 3.92

0.39 0.30 0.21 0.15 0.11 0.14

736.59 988.13 1468.36 2583.44 5717.25 13,871.21

2.65 2.78 3.26 4.25 5.25 6.80

1062.80 1677.52 3432.21 8764.47 24,149.36 81,703.92

2.66 2.78 3.09 3.69 4.02 3.88

0.29 0.21 0.14 0.11 0.14 0.23

1023.40 1542.84 2768.44 6129.98 14,468.47 25,098.95

2.56 2.58 2.61 2.67 2.83 3.41

518.49 637.05 824.75 1163.89 1931.33 4381.11

2.58 2.59 2.62 2.68 2.81 3.14

0.40 0.40 0.36 0.26 0.18 0.12

501.41 626.37 806.86 1117.14 1759.39 3471.83

2.55 2.56 2.58 2.61 2.65 2.76

403.05 471.40 567.39 711.70 951.67 1419.88

2.56 2.57 2.59 2.61 2.66 2.75

0.70 0.60 0.50 0.41 0.32 0.23

400.25 466.97 559.78 697.13 918.73 1324.60

0.5

0.1 0.2 0.3 0.4 0.5 0.6

2.62 2.67 2.79 3.21 4.13 5.10

766.44 1044.45 1615.72 3179.04 8102.23 23,945.97

2.63 2.68 2.80 3.05 3.60 3.94

0.39 0.30 0.21 0.15 0.11 0.14

750.91 1007.01 1495.18 2625.54 5788.69 14,014.82

2.67 2.80 3.28 4.27 5.26 6.81

1083.22 1708.22 3482.69 8845.65 24,300.82 82,014.36

2.68 2.80 3.09 3.71 4.05 3.91

0.29 0.21 0.15 0.11 0.14 0.23

1042.85 1570.79 2811.85 6201.47 14,608.79 25,353.75

2.58 2.60 2.63 2.69 2.86 3.43

528.76 649.63 840.93 1186.36 1966.55 4442.78

2.60 2.61 2.64 2.69 2.83 3.16

0.40 0.40 0.36 0.27 0.18 0.12

511.28 638.66 822.58 1138.49 1791.31 3526.25

2.57 2.58 2.60 2.63 2.67 2.78

411.05 480.75 578.61 725.70 970.16 1446.57

2.58 2.59 2.60 2.63 2.68 2.78

0.70 0.60 0.51 0.41 0.32 0.23

408.18 476.20 570.80 710.74 936.42 1349.25

0.75 0.1 0.2 0.3 0.4 0.5 0.6

2.64 2.69 2.82 3.23 4.15 5.11

781.04 1064.06 1644.78 3226.81 8181.60 24,105.17

2.65 2.70 2.81 3.08 3.62 3.97

0.39 0.30 0.22 0.15 0.11 0.14

765.12 1025.73 1521.78 2667.20 5859.36 14,156.78

2.69 2.83 3.31 4.29 5.28 6.83

1103.48 1738.63 3532.64 8926.30 24,451.56 82,323.89

2.70 2.82 3.11 3.73 4.07 3.94

0.29 0.21 0.15 0.11 0.14 0.23

1062.15 1598.49 2855.31 6272.18 14,747.15 25,604.92

2.60 2.62 2.65 2.71 2.88 3.46

538.96 662.12 856.99 1208.64 2001.45 4503.88

2.62 2.63 2.66 2.72 2.86 3.19

0.40 0.40 0.36 0.27 0.18 0.12

521.08 650.86 838.19 1159.69 1822.94 3580.12

2.59 2.60 2.62 2.65 2.69 2.80

419.00 490.02 589.74 739.58 988.50 1473.03

2.59 2.61 2.62 2.65 2.70 2.80

0.70 0.60 0.51 0.41 0.32 0.23

416.05 485.36 581.73 724.26 953.97 1373.68

Table 4 The variation of Rayleigh Number (Ra) with internal heat parameter (Q) at fixed values of Darcy–Brinkman number D˜ a = 0.4 and different values of anisotropy parameter (ζ ) for Q - source and γ - sink by linear/nonlinear analysis for case A, B, C, and D. D˜ a = 0.4 A

ζ

B

C

D

acL

RaL

acE

λ

RaE

acL

RaL

acE

λ

RaE

acL

RaL

acE

λ

RaE

acL

RaL

acE

λ

RaE

0.125 0.1 0.2 0.3 0.4 0.5 0.6

2.62 2.67 2.79 3.19 4.12 5.08

1447.75 1972.94 3053.85 6037.74 15,627.14 46,752.54

2.63 2.68 2.79 3.04 3.58 3.91

0.39 0.30 0.21 0.15 0.11 0.14

1418.78 1903.01 2827.76 4979.11 11,054.05 26,792.91

2.67 2.80 3.26 4.26 5.25 6.80

2046.51 3229.76 6624.51 17,101.25 47,517.77 161,807.39

2.68 2.80 3.11 3.69 4.01 3.87

0.29 0.21 0.14 0.11 0.14 0.23

1971.06 2971.60 5338.39 11,865.84 27,976.78 48,420.83

2.58 2.60 2.63 2.69 2.85 3.42

998.77 1226.98 1588.12 2240.49 3717.18 8463.02

2.60 2.61 2.64 2.70 2.83 3.15

0.40 0.40 0.36 0.26 0.18 0.12

965.84 1206.45 1553.94 2151.11 3387.77 6692.81

2.58 2.59 2.60 2.63 2.67 2.78

776.45 908.03 1092.79 1370.45 1832.02 2732.62

2.58 2.59 2.61 2.63 2.68 2.77

0.70 0.60 0.51 0.41 0.32 0.23

771.11 899.57 1078.23 1342.59 1769.09 2550.37

0.25 0.1 0.2 0.3 0.4

2.63 2.68 2.79 3.20

1455.09 1982.80 3068.50 6061.97

2.64 2.69 2.80 3.05

0.39 0.30 0.21 0.15

1425.92 1912.43 2841.16 50 0 0.21

2.68 2.80 3.27 4.26

2056.71 3245.10 6649.88 17,142.03

2.69 2.80 3.11 3.69

0.29 0.21 0.14 0.11

1980.78 2985.57 5360.13 11,901.95

2.59 2.61 2.64 2.70

1003.89 1233.26 1596.19 2251.70

2.61 2.62 2.64 2.70

0.40 0.40 0.36 0.26

970.76 1212.59 1561.78 2161.76

2.58 2.59 2.61 2.63

780.45 912.69 1098.38 1377.43

2.58 2.59 2.61 2.64

0.70 0.60 0.51 0.41

775.07 904.17 1083.73 1349.38

0.5 0.6 0.1 0.2 0.3 0.4 0.5 0.6

4.12 5.09 2.65 2.70 2.82 3.22 4.14 5.10

15,667.27 46,832.84 1484.31 2022.06 3126.77 6158.25 15,827.09 47,153.04

3.58 3.92 2.65 2.69 2.81 3.06 3.60 3.93

0.11 0.14 0.39 0.30 0.21 0.15 0.11 0.14

11,090.06 26,865.75 1440.18 1931.19 2867.87 5042.25 11,161.79 27,010.74

5.25 6.81 2.69 2.82 3.28 4.27 5.26 6.81

47,593.78 161,962.95 2077.03 3275.67 6700.42 17,223.36 47,745.51 162,273.74

4.01 3.88 2.70 2.82 3.13 3.71 4.03 3.89

0.14 0.23 0.29 0.21 0.14 0.11 0.14 0.23

28,048.23 48,550.75 20 0 0.14 3013.42 5403.43 11,973.86 28,190.45 48,809.09

2.86 3.42 2.60 2.62 2.65 2.71 2.87 3.43

3734.78 8494.00 1014.11 1245.78 1612.29 2274.04 3769.85 8555.74

2.84 3.16 2.62 2.63 2.65 2.71 2.85 3.17

0.18 0.12 0.40 0.40 0.36 0.27 0.18 0.12

3403.72 6720.10 980.59 1224.82 1577.43 2183.01 3435.50 6774.47

2.68 2.78 2.59 2.60 2.62 2.64 2.69 2.79

1841.25 2745.95 788.41 921.98 1109.54 1391.35 1859.64 2772.50

2.68 2.78 2.59 2.60 2.62 2.65 2.69 2.79

0.32 0.23 0.70 0.60 0.51 0.41 0.32 0.23

1777.92 2562.68 782.96 913.35 1094.69 1362.93 1795.52 2587.21

0.75 0.1 0.2 0.3 0.4 0.5 0.6

2.65 2.70 2.82 3.22 4.14 5.10

1484.31 2022.06 3126.77 6158.25 15,827.09 47,153.04

2.66 2.70 2.82 3.07 3.61 3.94

0.39 0.30 0.21 0.15 0.11 0.14

1454.37 1949.91 2894.46 5084.06 11,233.12 27,154.85

2.70 2.83 3.30 4.28 5.27 6.82

2097.27 3306.10 6750.70 17,304.43 47,896.87 162,584.07

2.71 2.83 3.14 3.72 4.04 3.91

0.29 0.21 0.14 0.11 0.14 0.23

2019.43 3041.14 5446.48 12,045.35 28,331.77 49,065.51

2.61 2.63 2.66 2.72 2.88 3.45

1024.30 1258.25 1628.32 2296.30 3804.76 8617.18

2.63 2.64 2.66 2.72 2.86 3.18

0.40 0.40 0.36 0.27 0.18 0.12

990.37 1237.00 1593.02 2204.19 3467.13 6828.55

2.60 2.61 2.63 2.65 2.70 2.81

796.34 931.25 1120.66 1405.22 1877.96 2798.95

2.60 2.61 2.63 2.66 2.71 2.80

0.70 0.60 0.51 0.41 0.32 0.23

790.82 922.50 1105.61 1376.42 1813.04 2611.63

0.5

ξ

tem is unstable and below the nonlinear plot, it is unstable. It is observed from these figures that as internal heat parameter increases, it tends to move the neutral stability curves upwards or in other sense the basic state tends to more stable and also, it occurs with the increase of anisotropy parameter. The Linear theory fails to capture the physics of fluid flow at a low value of wave number (a) but as wave number increases, two theories (linear and nonlinear) coincide. Figs. 7–10 displays the behavior of critical Rayleigh number with increasing values of internal heat parameter at different values of anisotropy parameter (ζ = 0.25, 0.75) and Darcy–Brinkman number (D˜ a = 0, 0.2) for linear and nonlinear analysis.

230

R. Nandal and A. Mahajan / Applied Mathematics and Computation 358 (2019) 216–231

From these figures, it appears that with an increase of internal heat parameter value of critical Rayleigh number increases and an increase in the value of anisotropy parameter further increases the value critical Rayleigh number, which results in the system more stable. It is noticed that increase in Darcy–Brinkman number highly raised the value of critical Rayleigh number and stabilize the system very much. It is further found that results of the linear and nonlinear analysis are very different which implies that the subcritical instability region exists. 6. Conclusion Here, the linear and nonlinear analysis is performed for a fluid saturated anisotropic porous medium with the internal heat generating function and extraction of heat through the boundary. From present work, it is concluded that heat extraction through boundary overcomes the effect of internal heat generation and enhances the stability of the system. Anisotropy parameter has to stabilize effect on the system because with an increase of its value of critical Rayleigh number increases. Results of the linear and nonlinear analysis are found different, so for all the heat generating function, the region of subcritical instability exists. It is found that increase in Darcy–Brinkman number delay the onset of convection and has a stabilizing effect on the system. Coupling parameter depends on the internal heat parameter and independent of the Darcy–Brinkman number and anisotropy parameter. References [1] A.A. Altawallbeh, B.S. Bhadauria, I. Hashim, Linear and nonlinear double-diffusive convection in a saturated anisotropic porous layer with Soret effect and internal heat source, Int. J. Heat Mass Trans. 59 (2013) 103–111. [2] K.A. Ames, S.S. Cobb, Penetrative convection in a porous medium with internal heat sources, Int. J. Engng. Sci 32 (1994) 95–105. [3] M. Berlengiero, K.A. Emanuel, J. von Hardenberg, A. Provenzale, E.A. Spiegel, Internally cooled convection: a fillip for Philip, Commun. Nonlinear Sci. Numer. Simul 17 (2012) 1998–2007. [4] B.S. Bhadauria, Double diffusive convection in a saturated anisotropic porous layer with internal heat source, Transp. Porous. Med. 92 (2012) 299–320. [5] A.N. Borujerdi, A.R. Noghrehabadi, D.A.S. Rees, Influence of Darcy number on the onset of convection in a porous layer with a uniform heat source, Int. J. Therm. Sci. 47 (2008) 1020–1025. [6] H.C. Brinkman, A calculation of the viscous force exerted by a flowing fluid on a dense swarm of particles, Appl. Sci. Res. 1 (1947) 27–34. [7] C. Canuto, M.Y. Hussaini, A. Quateroni, T. Zang, Spectral Methods in Fluid Dynamics, Springer, New York, 1998. [8] F. Capone, M. Gentile, A.A. Hill, Double-diffusive penetrative convection simulated via internal heating in an anisotropic porous layer with throughflow, Int. J. Heat Mass Trans. 54 (2011) 1622–1626. [9] F. Capone, M. Gentile, A.A. Hill, Penetrative convection in anisotropic porous media with variable permeability, Acta Mech 216 (2011) 49–58. [10] F. Capone, M. Gentile, A.A. Hill, Penetrative convection via internal heating in anisotropic porous media, Mech. Res. Commun. 37 (2010) 441–444. [11] M. Carr, B. Straughan, Penetrative convection in a fluid overlying a porous layer, Adv. Water Resour 26 (2003) 263–276. [12] G. Castinel, M. Combarnous, Natural convection in an anisotropic porous layer, Int. Chem. Engng. 17 (1977) 605–613. [13] S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability, Dover, New York, 1981. [14] G. Degan, C. Akowanou, L. Fagbemi, J. Zinsalo, Hydrodynamic anisotropy effects on radiation mixed convection interaction in a vertical porous channel, Appl. Math 7 (2016) 22–39. [15] J.E. Epherre, Critѐre d’apparition de la convection naturelle dans unee couche poreuse anisotrope, Rev. Thermique 168 (1975) 949–950. [16] S.N. Gaikwad, M. Dhanraj, The onset of double diffusive convection in a couple stress fluid saturated anisotropic porous layer with cross-diffusion effects, Int. J. Math. Arch 3 (2012) 2713–2727. [17] Y.H. Gandadharaiah, S.P. Suma, H. Nagarathnamma, Penetrative convection via internal heating in superposed fluid and anisotropic porous layers with throughflow, Int. J. Math. Comp. Sci. 3 (2014) 5–28. [18] D. Goluskin, E.A. Spiegel, Convection driven by internal heating, Phys. Lett. A 377 (2012) 83–92. [19] D. Goluskin, E.P. van der Poel, Penetrative internally heated convection in two and three dimensions, J. Fluid Mech 791 (2016) 1–13. [20] P. Guo, Dependency of tortuosity and permeability of porous media on directional distribution of pore voids, Trans. Porous Med. 95 (2012) 285–303. [21] A.J. Harfash, Stability analysis of penetrative convection in anisotropic porous media with variable permeability, J. Non Equilib. Thermodyn. 39 (2014) 123–133. [22] T. Hattori, J.C. Patterson, C. Lei, On the stability of transient penetrative convection driven internal heating coupled with a thermal boundary condition, Int. Comm. Heat Mass Transf. 64 (2015) 29–33. [23] I.H. Herron, Onset of convection in a porous medium with internal heat source and variable gravity, Int. J. Eng. Sci. 39 (2001) 201–208. [24] E. Kaminski, A. Chenet, C. Jaupart, V. Courtillot, Rise of volcanic plumes to the stratosphere aided by penetrative convection above large lava flows, Earth Planetary Sci 301 (2011) 171–178. [25] T. Karmakar, G.P.R. Sekhar, Effects of anisotropic permeability on fluid flow through composite channel, J. Eng. Math. 100 (2016) 33–51. [26] A.V. Kuznetsov, D.A. Nield, The effect of spatially nonuniform internal heating on the onset of convection in a horizontal fluid layer, J. Heat Trans. 138 (2016) 1–9. [27] A. Mahajan, M.K. Sharma, Penetrative convection in magnetic nanofluids via internal heating, Phys. Fluids 1–14 (2017) 29. [28] A. Mahajan, R. Nandal, Anisotropic porous penetrative convection for a local thermal non-equilibrium model with Brinkman effects, Int. J. Heat Mass Trans. 115 (2017) 235–250. [29] A. Mahidjiba, L. Robillard, P. Vasseur, M. Mamou, Onset of convection in an anisotropic porous layer of finite lateral extent, Int. Comm. Heat Mass Trans. 27 (20 0 0) 333–342. [30] M.S. Malashetty, P. Kollur, The onset of double diffusive convection in a couple stress fluid saturated anisotropic porous layer, Trans. Porous Med 86 (2011) 435–459. [31] P.C. Matthews, A model for the onset of penetrative convection, J. Fluid Mech. 188 (1988) 571–583. [32] P.C. Matthews, A model for the onset of penetrative convection, J. Fluid Mech. 188 (1988) 571–583. [33] R. Nandal, A. Mahajan, Linear and nonlinear stability analysis of a Horton-Rogers-Lapwood problem with an internal heat source and Brinkman effects, Transp. Porous Med. 117 (2017) 261–280. [34] W.F. Orr, Stability or instability of the steady motions of a perfect liquid, in: Proceedings of the Royal Irish Academy, 1907 27, 9-68 and 69-138. [35] C. Parthiban, P.R. Patil, Thermal instability in an anisotropic porous medium with internal heat source and inclined temperature gradient, Int. Comm. Heat Mass Trans. 24 (1997) 1049–1058. [36] M. Prud’homme, S. Jasmin, Inverse solution for a biochemical heat source in a porous medium in the presence of natural convection, Chem. Engg. Sci. 61 (2006) 1667–1675. [37] Y. Qin, P.N. Kaloni, Convective instabilities in anisotropic porous media, Stud. Appl. Math. 91 (1994) 189–204. [38] O. Reynolds, On the dynamical theory of incompressible viscous fluids and the determination of the criterion, Philosoph. Trans. R. Soc. Lond. 186 (1895) 123–164.

R. Nandal and A. Mahajan / Applied Mathematics and Computation 358 (2019) 216–231 [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52]

231

P.A. Rice, D.J. Fontugne, R.G. Latini, A.J. Barduhn, Anisotropic permeability in porous media, Ind. Eng. Chem. 62 (1970) 23–31. P.H. Roberts, Convection in horizontal layers with internal heat generation, Theory J. Fluid Mech. 30 (1967) 33–49. S. Safi, S. Benissaad, Heat and mass transfer in anisotropic porous media, Adv. Theor. Appl. Mech 5 (2012) 11–22. I.S. Shivakumara, M. Dhananjaya, Penetrative Brinkman convection in an anisotropic porous layer saturated by a nanofluid, Ain Shams Eng. J. 6 (2015) 703–713. I.S. Shivakumara, S.P. Suma, R. Indira, Y.H. Gandadharaiah, Effect of internal heat generation on the onset of Marangoni convection in a fluid layer overlying a layer of an anisotropic porous medium, Transp. Porous Med. 92 (2012) 727–743. P.K. Srimani, H.R. Sudhakar, Linear and nonlinear penetrative convection in a porous layer, Indian J. Pure Appl. Math 27 (1996) 499–518. L. Storesletten, Natural convection in a horizontal porous layer with anisotropic thermal diffusivity, Trans. Porous Med. 12 (1993) 19–29. B. Straughan, D.W. Walker, Anisotropic porous penetrative convection, Proc. R. Soc. Lond. A 452 (1996) 97–115. B. Straughan, Continuous dependence on the heat source and non-linear stability in penetrative convection, Int. J. Non Linear Mech. 26 (1991) 221–231. B. Straughan, Continuous dependence on the heat source and non-linear stability for convection with internal heat generation, Math. Methods Appl. Sci. 13 (1990) 373–383. B. Straughan, Resonant penetrative convection with an internal heat source/sink, Acta Appl. Math. 132 (2014) 561–581. B. Straughan, The Energy Method, Stability and Nonlinear Convection, Springer-Verlag, New York, 2004. P.A. Tyvand, L. Storesletten, Onset of convection in an anisotropic porous medium with oblique principal axes, J. Fluid Mech. 226 (1991) 371–382. A.P. Van den Berg, D.A. Yuen, G.L. Beebe, M.D. Christiansen, The dynamical impact of electronic thermal conductivity on deep mantle convection of exosolar planets, Phys. Earth Planet Inter. 178 (2010) 136–154.