# Heat Transfer and Flows of Thermal Convection in a Fluid-Saturated Rotating Porous Medium

Authors: Jianhong Kang; Tongqiang Xia; Yingke Liu
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Thermal convection at the steady state for high Rayleigh number in a rotating porous half space is investigated. Taking into account the effect of rotation, Darcy equation is extended to incorporate the Coriolis force term in a rotating reference frame. The velocity and temperature fields of thermal convection are obtained by using the homotopy analysis method. The influences of Taylor number and Rayleigh number on the Nusselt number, velocity profile, and temperature distribution are discussed in detail. It is found that the Nusselt number decreases rapidly with the increase of Taylor number but tends to have an asymptotic value. Besides, the rotation can give rise to downward flow in contrast with the upward thermal convection.


---

## Body

## 1. Introduction

Thermal convection in porous media has been given a great deal of attention by researchers. Its importance has relevance to a wide range of phenomena in geophysical, astrophysical, and engineering applications [1–3]. Since the pioneering work of Horton and Rogers [4], thermal convection in porous media has been investigated extensively. Lapwood [5] and Katto and Masuoka [6] examined thermal convective instability in a fluid-saturated horizontal porous layer. Wang [7] investigated thermal convective instability in a fluid-saturated rectangular box. Hayat et al. [8] and Yin et al. [9] considered thermal convective stability of a viscoelastic fluid in a porous layer. Malashetty et al. [10] considered thermal convective instability in an anisotropic porous layer. Nield and Kuznetsov [11] studied thermal convective instability in a nanofluid-saturated porous layer. All of these references focused on determining the criterion for the onset of thermal convection in terms of the critical Rayleigh number. Only if the Rayleigh number exceeds the critical value, thermal convection can set in; otherwise, heat transfer is dominated by conduction [12]. At the onset of thermal convection (i.e., the low-Rayleigh-number convection), the deviation of temperature and velocity from the purely conductive basic state is considered to be infinitesimal [13]. Consequently, the nonlinear terms in the governing equations are higher-order infinitesimals compared with other linear terms [14]. For the first-order approximation, the nonlinear terms can be neglected such that the negligence of nonlinearity of the problem is justified [15]. However, with the increase of Rayleigh number, the amplitudes of thermal convection grow rapidly such that the nonlinearity cannot be neglected [16, 17]. Due to the difficulty arising from the nonlinearity, only few works with numerical methods have been done [18, 19]. Till now, the study on high-Rayleigh-number convection in porous media is far less fruitful.A particularly interesting variation of thermal convection is the case where the sample is rotated about a vertical axis with a uniform angular speed. In this system, the Coriolis and centrifugal forces besides gravity should be taken into account [20]. The Coriolis effect on thermal convection in different rotating systems has been studied extensively in the literature. Chandrasekhar [21] analyzed the thermal convective stability in a rotating horizontal layer of pure fluid, which was extended to the case of porous layer by Vadasz [22] and then the case of viscoelastic fluid by Kang et al. [23]. Malashetty et al. [24, 25] also considered thermal stability in a rotating anisotropic porous layer and in a rotating porous layer by using a thermal nonequilibrium model. Jones et al. [26] studied linear stability of thermal convection in a rapidly rotating sphere and developed the asymptotic theory for the sphere, which was later extended to the case of spherical shells by Dormy et al. [27]. Kang et al. [28] conducted linear stability analysis of viscoelastic fluids in a rotating porous cylindrical annulus. Besides, Busse [29] studied the onset of convection in a narrow cylindrical annulus heated from below and rotating about its vertical axis of symmetry. Ponty et al. [30] studied the onset of thermal convection in a rotating horizontal layer where the rotating axis was inclined at an angle to the vertical. Straughan [31] performed a sharp nonlinear stability threshold in rotating porous convection. These works are based on the assumption that not far from the axis of rotation gravity buoyancy is dominant and the centrifugal buoyancy can be neglected [22, 29]. On the other hand, the centrifugal buoyancy-driven thermal convection was also investigated where gravity buoyancy was assumed to be negligible compared with the centrifugal buoyancy [32, 33]. Assuming a small top aspect ratio of a narrow porous layer, Vadasz [34, 35] originally investigated thermal convection induced by centrifugal force in a rotating porous layer. His work was later extended by many researchers to consider other effects on the centrifugally driven thermal convection in porous media. For example, the anisotropic effects of porous media were studied by Govender [36] and Saravanan and Brindha [37]; the effects of thermal modulation and rotation modulation were considered by Suthar and Bhadauria [38]. Besides, Alloui and Vasseur [39] studied convective heat transfer with centrifugal force field for a horizontal porous annulus of infinite extent. Kang et al. [40] investigated thermal convection driven by centrifugal buoyancy in a rotating porous cylindrical annulus with conical end boundaries. However, to the best of the authors’ knowledge, no paper has investigated the Coriolis effect on the heat transfer rate for the high-Rayleigh-number thermal convection in porous media, which is a coupled nonlinear problem.In general, it is difficult to solve a nonlinear problem analytically and any analytical solution for nonlinear partial differential equations is always exciting. Several analytical techniques including perturbation, artificial small parameter, theδ-expansion method, variational iteration, and Adomian decomposition methods have been developed to solve nonlinear problems [41, 42], among which the perturbation method and Adomian decomposition method are most widely applied. The perturbation method transforms a nonlinear problem into a family of infinite linear problems, and the superimposed solutions of linear problems are used to approximate the solution for nonlinear problem [43]. However, perturbation techniques are essentially based on small perturbation quantities, and the so-called “small parameter assumption” of perturbation techniques greatly restricts their applications [44]. Furthermore, the perturbation technique may become invalid for strongly nonlinear problems [45]. In contrast, the Adomian decomposition method is an effective method to solve strongly nonlinear problem and does not depend on any small parameter assumption [46]. The Adomian decomposition method is based on the assumption that the differential operator of nonlinear problems can be divided into linear and nonlinear parts; then the solution for nonlinear problem may be expressed as a polynomial series [47]. The polynomial series obtained by Adomian decomposition method converges quickly, but its convergence radius is relatively small [48]. More recently, Liao [49] developed a new kind of analytic method for nonlinear problems, called the homotopy analysis method (HAM), which is based on the homotopy, a basic concept of topology. Different from perturbation techniques, the homotopy analysis method is valid even for nonlinear problems whose governing equations and/or boundary conditions do not contain any small parameters at all [44]. The homotopy analysis method also provides us with great freedom to select proper base functions to approximate solutions of nonlinear problems [50]. Furthermore, the homotopy analysis method provides us with a family of solution series and a simple way to adjust and control the convergence region and rate of approximation series [48]. Serving as a powerful tool to deal with nonlinear equations, the homotopy analysis method has applied to many nonlinear problems in science and engineering [51, 52]. Abbasbandy [53] used the homotopy analysis method to solve a generalized Hirota-Satsuma coupled KdV equation. Hayat et al. [54] investigated mixed convection flow of a micropolar fluid over a nonlinearly stretching sheet by using the homotopy analysis method. Mehmood et al. [55] presented complete analytic solution to the unsteady heat transfer flow of an incompressible viscous fluid over a permeable plane wall with the homotopy analysis method.In this work, we study high-Rayleigh-number steady state thermal convection in a rotating porous half space. Our focus is how an external constraint of rotation can affect the flow features and heat transfer in a porous half space. The present analysis combines the effects of high Rayleigh number and rotation together. We solve the coupled nonlinear equations by using the homotopy analysis method (HAM). The other parts of this paper are organized as follows. In Section2, we give a mathematical formulation of the considered problem. In Section 3, the series solutions for velocity field and temperature distribution are presented using the HAM. Section 4 is devoted to main results and discussion on the convective flow and heat transfer in a rotating porous medium. Finally, this paper closes with a conclusion in Section 5.


## 2. Mathematical Formulation

Consider a porous half space saturated with a viscous fluid with the densityρ and viscosity μ, as shown in Figure 1. The porous medium with permeability K and porosity ϕ is heated from below with constant temperature T0 and cooled from above at infinity with constant temperature T∞. Meanwhile, the porous medium rotates around z-axis with uniform angular velocity Ω. Here, a rotating cylindrical coordinate frame of reference is chosen to link together with the porous matrix. Not far (at distances d≪g/Ω2) from the axis of rotation, one can justify that the gravity buoyancy is dominant and centrifugal buoyancy can be neglected [22].Figure 1
Schematic diagram of the physical problem.In modeling flows in a fluid-saturated porous medium, the simplest but most widely used model in engineering is the well-known Darcy model which states proportionality between flow rate and the pressure gradient and is generally valid in a densely packed porous medium of sufficiently large thickness [1, 2]. For a loosely packed porous media, the Darcy-Brinkman model is well accepted [56], which added a viscous term so that the no-slip condition on solid surfaces can be applied [57, 58]. When the thermal convective instability in porous media is considered, a modified Darcy model which extended the Darcy model to include a time derivative term is usually adopted to study the oscillatory convection as well as stationary convection [22, 40, 59]. Besides, by analogy with the constitutive equation of Oldroyd-B fluid, some researchers proposed a Darcy-Oldroyd model for modeling flows of viscoelastic fluids in porous media, which introduces a relaxation parameter characterizing the elasticity of viscoelastic fluid and a retardation parameter characterizing the viscous damping of viscoelastic fluid [13, 60, 61].In this work, the Darcy model is extended with the addition of the Coriolis term in rotating reference frame and is applied to characterize the flow in porous media [1]. Assuming a state of local thermal equilibrium between the fluid and solid phase, the steady state momentum, energy, and continuity equations for this system are given as [17, 18](1)μKV+2ρΩϕk×V=-∇P-ρgk,V·∇T=κ∇2T,∇·V=0in which V,P, and T are the Darcian velocity, pressure, and temperature. g is the gravitational acceleration, k is a unit vector in the vertical direction, and κ is the effective thermal diffusivity, respectively.Boussinesq (1903) and Oberbeck (1879) had recognized that when temperature variations are small, the buoyancy of fluid is significant and they assumed an equation of state in which the density is a linear function of temperature and independent of pressure [15]; that is,(2)ρ=ρ01-αT-T0,where ρ0 is the density at temperature T0 and α is the volumetric expansion coefficient of the fluid. The assumption stated in (2) is known as Boussinesq approximation and is widely adopted in the study on thermal convection [5–19, 21–28].Chooseκ/d, μκ/K, and T0-T∞ to scale the velocity, pressure, and temperature respectively; we can obtain from (1)-(2) the dimensionless governing equations in the following form:(3)V+Ta1/2k×V=-∇P+Raθk,(4)V·∇θ=∇2θ,(5)∇·V=0,where we have used symbol θ=(T-T∞)/(T0-T∞) to denote dimensionless temperature and kept the same notation for all other dimensionless variables for brevity. The dimensionless parameters introduced in (3) are Ta, the Taylor number, and Ra, the Rayleigh number, defined, respectively, as(6)Ta=2ρ0ΩKμϕ2,Ra=ρ0gαT0-T∞Kdμκ.Previous studies have indicated that two-dimensional flow structures of thermal convection are possible in a fluid-saturated porous layer [16–19]. Due to the existence of rotation axis, the two-dimensional assumption of azimuthal symmetry is made in the present paper. Dropping the azimuthal dependence of all variables in (3)–(5), a stream function ψ is introduced such that (u,w)=((1/r)(∂ψ/∂z),-(1/r)(∂ψ/∂r)), where u and w are the components of V in the r- and z-directions, respectively. Then, (3)–(5) are reduced to the following form:(7)1+Ta∂2ψ∂z2+∂2ψ∂r2-1r∂ψ∂r+rRa∂θ∂r=0,1r∂ψ∂z∂θ∂r-1r∂ψ∂r∂θ∂z=∂2θ∂z2+∂2θ∂r2+1r∂θ∂r.The boundary conditions for (7) can be formulated as (i) θ=1 and ψ=0, at z=0; (ii) θ=0 and ∂ψ/∂z=0, at z→∞.By introducing the stream function, the considered problem is greatly simplified so that it can be solved analytically by the homotopy analysis method (HAM).


## 3. Analytic Solutions by the Homotopy Analysis Method

In this section, we intend to obtain series solutions for the coupled nonlinear equations (7) by virtue of the homotopy analysis method (HAM) [50, 62]. For this purpose, the following transformations are introduced: ψ=f(r,η) and θ=φ(r,η) where η=z/r. Then (7) are readily transformed into the following equations:(8)∂2f∂r2+1+Ta+η2r2∂2f∂η2-2ηr∂2f∂η∂r+3ηr2∂f∂η-1r∂f∂r+Rar∂φ∂r-η∂φ∂η=0,∂2φ∂r2+1+η2r2∂2φ∂η2-2ηr∂2φ∂η∂r+ηr2∂φ∂η+1r∂φ∂r+1r2∂f∂r∂φ∂η-∂f∂η∂φ∂r=0.The corresponding boundary conditions are (i)φ=1 and f=0, at η=0; (ii) φ=0 and ∂f/∂η=0, at η→∞.


### 3.1. The Zero-Order Deformation Equations

Following the basic idea of HAM, we assume thatf(η,r) and g(η,r) can be expressed by a set of base functions {riηje-kγη∣k≥0,j≥0,i≥0}:(9)fr,η=∑k=0∞∑j=0∞∑i=0∞ai,jkriηje-kγη,φr,η=∑k=0∞∑j=0∞∑i=0∞bi,jkriηje-kγη,where ai,jk and bi,jk are coefficients and γ is a spatial-scale parameter.According to the solution expressions (9) together with corresponding boundary conditions, it is straightforward to choose the following auxiliary linear operators and initial guess solutions of f(η,r) and g(η,r):(10)LfΦr,η;q=∂2Φ∂η2-γ2Φ,(11)LφΘr,η;q=∂2Θ∂η2+γ∂Θ∂η,(12)f0r,η=1-e-γη,φ0r,η=e-γη.Now, we construct the so-called zero-deformation equations as follows:(13)1-qLfΦr,η;q-f0r,η=qħHr,ηNfΦr,η;q,q∈0,1,1-qLφΘr,η;q-φ0r,η=qħHr,ηNφΘr,η;q,q∈0,1subject to the boundary conditions(14)Φr,0;q=0,Θr,0;q=1,∂Φr,η;q∂η∂Φr,η;q∂ηη→∞=0,Θr,η;qr,η;qη→∞=0,where q is the embedding parameter, ħ the nonzero auxiliary parameter, and H(r,η) the nonzero auxiliary function in HAM. Since the auxiliary function would not have effect on the final solution, it is assumed that H(r,η)=1 in the further discussion. Nf and Nφ are nonlinear operators defined as(15)NfΦr,η;q=1+Ta+η2∂2Φ∂η2-2rη∂2Φ∂η∂r+r2∂2Φ∂r2+3η∂Φ∂η-r∂Φ∂r+r2Rar∂Θ∂r-η∂Θ∂η,NφΘr,η;q=1+η2∂2Θ∂η2-2rη∂2Θ∂η∂r+r2∂2Θ∂r2+η∂Θ∂η+r∂Θ∂r+∂Φ∂r∂Θ∂η-∂Φ∂η∂Θ∂r.Obviously, whenq=0, the solutions of the zero-order deformation equations (13) reduce to the initial guess solutions(16)Φr,η;0=f0r,η,Θr,η;0=φ0r,η;when q=1, (13)-(14) are exactly the same as (9)–(11); that is,(17)Φr,η;1=fr,η,Θr,η;1=φr,η.It can be observed that asq increases from 0 to 1, Φ(r,η;q) and Θ(r,η;q) vary continuously from the initial guess solutions f0(r,η) and φ0(r,η) to the exact solutions f(r,η) and φ(r,η) for the considered problem. This kind of continuous variation is called deformation in topology.Now, we expandΦ(r,η;q) and Θ(r,η;q) in a power series of the embedding parameter q as follows:(18)Φr,η;q=f0r,η+∑k=0∞‍fmr,ηqm,Θr,η;q=φ0r,η+∑k=0∞‍φmr,ηqm,where (19)fm(r,η)=1m!∂mΦ∂qm∂mΦ∂qmq=0,φmr,η=1m!∂mΘ∂qm∂mΘ∂qmq=0.It is noteworthy that (10)–(13) introduce two nonzero parameters ħ and γ so that the series expressions (18) depend on the two parameters. The auxiliary parameters ħ and γ play important roles in the HAM, which are used to guarantee the convergence of the final series solutions. In general, ħ and γ are chosen so that series (18) are convergent at q=1. Therefore, the analytic solutions of the problem by HAM are given as(20)fr,η=f0r,η+∑m=0∞‍fmr,η,φr,η=φ0r,η+∑m=0∞‍φmr,η.


### 3.2. The High-Order Deformation Equations

In (20), the functions fm and φm (m=1,2,3,…) remain unknown. Differentiating the zero-order deformation equations (13) m times with respect to q and setting q=0, the m-order deformation equations can be obtained as(21)Lffmr,η-χmfm-1r,η=ħRmfm-1,fm-2,…,f0,Lφfmr,η-χmφm-1r,η=ħQmφm-1,φm-2,…,φ0subject to boundaries conditions(22)fmr,0=0,φmr,0=0,∂fm(r,η)∂η∂fm(r,η)∂ηη→∞=0,φmr,ηr,ηη→∞=0,where (23)χm=0,m=11,m>1,(24)
R
m
f
m
-
1
,
…
,
f
0
=
1
m
-
1
!
∂
m
-
1
N
f
Φ
r
,
η
;
q
∂
q
m
∂
m
-
1
N
f
Φ
r
,
η
;
q
∂
q
m
q
=
0
=
1
+
Ta
+
η
2
∂
2
f
m
-
1
∂
η
2
-
2
r
η
∂
2
f
m
-
1
∂
η
∂
r
+
r
2
∂
2
f
m
-
1
∂
r
2
+
3
η
∂
f
m
-
1
∂
η
-
r
∂
f
m
-
1
∂
r
+
r
2
Ra
r
∂
φ
m
-
1
∂
r
-
η
∂
φ
m
-
1
∂
η
,
Q
m
φ
m
-
1
,
…
,
φ
0
=
1
m
-
1
!
∂
m
-
1
N
φ
Θ
r
,
η
;
q
∂
q
m
∂
m
-
1
N
φ
Θ
r
,
η
;
q
∂
q
m
q
=
0
=
1
+
η
2
∂
2
φ
m
-
1
∂
η
2
-
2
r
η
∂
2
φ
m
-
1
∂
η
∂
r
+
r
2
∂
2
φ
m
-
1
∂
r
2
+
η
∂
φ
m
-
1
∂
η
+
r
∂
φ
m
-
1
∂
r
+
∑
i
=
0
m
-
1
‍
∂
f
i
∂
r
∂
φ
m
-
i
-
1
∂
η
-
∂
f
i
∂
η
∂
φ
m
-
i
-
1
∂
r
.Letfm* and φm* denote particular solutions of (21); then we can express the general solutions as follows:(25)fm=C1e-γη+C2eγη+fm*,φm=C3+C4e-γη+φm*,where the coefficients C1,C2,C3, and C4 can be determined by the boundary conditions (22). By virtue of symbolic computation software, one can compute fm and φm from (21)–(24) order by order.


## 3.1. The Zero-Order Deformation Equations

Following the basic idea of HAM, we assume thatf(η,r) and g(η,r) can be expressed by a set of base functions {riηje-kγη∣k≥0,j≥0,i≥0}:(9)fr,η=∑k=0∞∑j=0∞∑i=0∞ai,jkriηje-kγη,φr,η=∑k=0∞∑j=0∞∑i=0∞bi,jkriηje-kγη,where ai,jk and bi,jk are coefficients and γ is a spatial-scale parameter.According to the solution expressions (9) together with corresponding boundary conditions, it is straightforward to choose the following auxiliary linear operators and initial guess solutions of f(η,r) and g(η,r):(10)LfΦr,η;q=∂2Φ∂η2-γ2Φ,(11)LφΘr,η;q=∂2Θ∂η2+γ∂Θ∂η,(12)f0r,η=1-e-γη,φ0r,η=e-γη.Now, we construct the so-called zero-deformation equations as follows:(13)1-qLfΦr,η;q-f0r,η=qħHr,ηNfΦr,η;q,q∈0,1,1-qLφΘr,η;q-φ0r,η=qħHr,ηNφΘr,η;q,q∈0,1subject to the boundary conditions(14)Φr,0;q=0,Θr,0;q=1,∂Φr,η;q∂η∂Φr,η;q∂ηη→∞=0,Θr,η;qr,η;qη→∞=0,where q is the embedding parameter, ħ the nonzero auxiliary parameter, and H(r,η) the nonzero auxiliary function in HAM. Since the auxiliary function would not have effect on the final solution, it is assumed that H(r,η)=1 in the further discussion. Nf and Nφ are nonlinear operators defined as(15)NfΦr,η;q=1+Ta+η2∂2Φ∂η2-2rη∂2Φ∂η∂r+r2∂2Φ∂r2+3η∂Φ∂η-r∂Φ∂r+r2Rar∂Θ∂r-η∂Θ∂η,NφΘr,η;q=1+η2∂2Θ∂η2-2rη∂2Θ∂η∂r+r2∂2Θ∂r2+η∂Θ∂η+r∂Θ∂r+∂Φ∂r∂Θ∂η-∂Φ∂η∂Θ∂r.Obviously, whenq=0, the solutions of the zero-order deformation equations (13) reduce to the initial guess solutions(16)Φr,η;0=f0r,η,Θr,η;0=φ0r,η;when q=1, (13)-(14) are exactly the same as (9)–(11); that is,(17)Φr,η;1=fr,η,Θr,η;1=φr,η.It can be observed that asq increases from 0 to 1, Φ(r,η;q) and Θ(r,η;q) vary continuously from the initial guess solutions f0(r,η) and φ0(r,η) to the exact solutions f(r,η) and φ(r,η) for the considered problem. This kind of continuous variation is called deformation in topology.Now, we expandΦ(r,η;q) and Θ(r,η;q) in a power series of the embedding parameter q as follows:(18)Φr,η;q=f0r,η+∑k=0∞‍fmr,ηqm,Θr,η;q=φ0r,η+∑k=0∞‍φmr,ηqm,where (19)fm(r,η)=1m!∂mΦ∂qm∂mΦ∂qmq=0,φmr,η=1m!∂mΘ∂qm∂mΘ∂qmq=0.It is noteworthy that (10)–(13) introduce two nonzero parameters ħ and γ so that the series expressions (18) depend on the two parameters. The auxiliary parameters ħ and γ play important roles in the HAM, which are used to guarantee the convergence of the final series solutions. In general, ħ and γ are chosen so that series (18) are convergent at q=1. Therefore, the analytic solutions of the problem by HAM are given as(20)fr,η=f0r,η+∑m=0∞‍fmr,η,φr,η=φ0r,η+∑m=0∞‍φmr,η.


## 3.2. The High-Order Deformation Equations

In (20), the functions fm and φm (m=1,2,3,…) remain unknown. Differentiating the zero-order deformation equations (13) m times with respect to q and setting q=0, the m-order deformation equations can be obtained as(21)Lffmr,η-χmfm-1r,η=ħRmfm-1,fm-2,…,f0,Lφfmr,η-χmφm-1r,η=ħQmφm-1,φm-2,…,φ0subject to boundaries conditions(22)fmr,0=0,φmr,0=0,∂fm(r,η)∂η∂fm(r,η)∂ηη→∞=0,φmr,ηr,ηη→∞=0,where (23)χm=0,m=11,m>1,(24)
R
m
f
m
-
1
,
…
,
f
0
=
1
m
-
1
!
∂
m
-
1
N
f
Φ
r
,
η
;
q
∂
q
m
∂
m
-
1
N
f
Φ
r
,
η
;
q
∂
q
m
q
=
0
=
1
+
Ta
+
η
2
∂
2
f
m
-
1
∂
η
2
-
2
r
η
∂
2
f
m
-
1
∂
η
∂
r
+
r
2
∂
2
f
m
-
1
∂
r
2
+
3
η
∂
f
m
-
1
∂
η
-
r
∂
f
m
-
1
∂
r
+
r
2
Ra
r
∂
φ
m
-
1
∂
r
-
η
∂
φ
m
-
1
∂
η
,
Q
m
φ
m
-
1
,
…
,
φ
0
=
1
m
-
1
!
∂
m
-
1
N
φ
Θ
r
,
η
;
q
∂
q
m
∂
m
-
1
N
φ
Θ
r
,
η
;
q
∂
q
m
q
=
0
=
1
+
η
2
∂
2
φ
m
-
1
∂
η
2
-
2
r
η
∂
2
φ
m
-
1
∂
η
∂
r
+
r
2
∂
2
φ
m
-
1
∂
r
2
+
η
∂
φ
m
-
1
∂
η
+
r
∂
φ
m
-
1
∂
r
+
∑
i
=
0
m
-
1
‍
∂
f
i
∂
r
∂
φ
m
-
i
-
1
∂
η
-
∂
f
i
∂
η
∂
φ
m
-
i
-
1
∂
r
.Letfm* and φm* denote particular solutions of (21); then we can express the general solutions as follows:(25)fm=C1e-γη+C2eγη+fm*,φm=C3+C4e-γη+φm*,where the coefficients C1,C2,C3, and C4 can be determined by the boundary conditions (22). By virtue of symbolic computation software, one can compute fm and φm from (21)–(24) order by order.


## 4. Results and Discussion

An important dimensionless quantity in engineering applications is the area-averaged Nusselt number [63], which measures the heat transfer rate and can be defined as(26)Nu=Qz=0λΔT/dA,where 〈·〉A denotes area-averaged integration over a circular area A, Q|z=0 denotes the heat flux at the bottom of porous media, and λ is the thermal conductivity. The denominator of (26) can be taken as a heat flux by conduction in a fluid layer with the depth d and temperature difference ΔT at the top and bottom, and in the definition of Nusselt number it serves as a reference of heat flux to measure the convective heat transfer rate of thermal convection. The Nusselt number in terms of dimensionless variables can be rewritten as (27)Nu=∂θ∂z∂θ∂zz=0A=2∫01∂φ∂η∂φ∂ηη=0dr.


### 4.1. The Discussion on the Parametersħ and γ

The discussed problem is governed by coupling nonlinear partial differential equations (3)–(5), which cannot be analytically solved by traditional methods. In this work, analytical solutions are presented by virtue of the homotopy analysis method (HAM) based on the homotopy, a basic concept of topology. Unlike previous analytic techniques, the HAM provides us with a family of solution series by introducing the nonzero auxiliary parameter ħ. The parameter ħ plays an important role in the HAM. However, it cannot be chosen arbitrary because the convergence of solution series is dependent on its value.In order to clarify how to determine the valid value or interval of parameterħ, the so-called ħ-curves in the HAM are used, in which some related physical quantities versus parameter ħ are plotted. According to Liao [50, 62] who first proposed the HAM, if there exists a nearly horizontal line segment in the ħ-curve, the horizontal line segment corresponds to the valid interval of ħ, and the solution series are convergent for each ħ in the valid interval.In the present study, the curve Nu versusħ serves as the ħ-curve. Our strategy for determining the valid value of ħ proceeds as follows. Firstly, for a set of given parameters Ra,Ta, and γ, we plot the Nu-ħ curve. Secondly, we check whether there is a nearly horizontal line segment in the Nu-ħ curve. If there is none, the choice of γ is not appropriate; we change the value of γ and replot the Nu-ħ curve until the horizontal line segment appears. Finally, the valid interval of ħ and the value of γ are determined accordingly. Repeating this procedure for different Ra and Ta, we can determine all the valid values of ħ and γ used in our computation.Using the procedure mentioned above, the Nu-ħ curve with Ra=1000 and Ta=4 is shown in Figure 2. It can be seen that, for γ=9, there exist nearly horizontal line segments, and the valid interval of ħ is about [-0.3,-0.5]. For each value of ħ in the interval [-0.3,-0.5], the solution series (20) and (27) of the problem are expected to be convergent. However, it is found that, for other values of γ, no such horizontal line segment exists, indicating that, in those cases, the solution series are divergent for any ħ. Figure 3 shows the effect of valid ħ on the temperature distribution, where three values of ħ in the valid interval are taken. It can be seen that the difference among three curves is very small, which implies that if ħ is chosen from the valid interval its effect on the solution series is negligible.Figure 2
The display of Nu-ħ curve which is used to determine the valid interval of ħ.Figure 3
The effect of validħ on the temperature distribution.It is noteworthy that, in practical applications, the infinite solution series are truncated when the subscript satisfiesm>N (N is a positive integer). The accuracy of the truncation scheme is checked by choosing sufficiently large N such that the solutions do not change significantly when N is replaced by N+2. The Nu-ħ curves with Ra=1000, Ta=4, and γ=9 for different truncated order N are shown in Figure 4. It can be seen that the 11th-order truncated approximation has already been satisfactory.Figure 4
The Nu-ħ curves for different truncated order N with Ra=1000, Ta=4, and γ=9.In Table1, the values of Nusselt number with Ra=1000 and Ta=4 are calculated for different valid parameter ħ and truncated order N. From this table, it can be seen that the convergence rate of solution series for different valid ħ is different. Therefore, the parameter ħ can not only control whether the solution series are convergent but also adjust the convergence rate of solution series. Furthermore, it is easily obtained that the relative error between 11th-order truncated approximation and 9th-order truncated approximation is smaller than 1, which is accurate enough to predict the heat transfer rate.Table 1
The values of Nusselt number Nu for different truncated orderN and valid parameter ℏ of solution series with Ra=1000 and Ta=4.
Truncated orderN
The values of Nusselt number Nu
ℏ
=
-
0.30
ℏ
=
-
0.35
ℏ
=
-
0.40
ℏ
=
-
0.45
7
3.73733
3.38959
3.27674
3.25194
8
3.46868
3.29490
3.27273
3.27276
9
3.37163
3.27839
3.26864
3.26225
10
3.34877
3.26358
3.26591
3.26846
11
3.27582
3.26235
3.26534
3.26538The valid values of parametersħ and γ used in the following computation are partly listed in Table 2. It can be seen that the value of γ tends to increase with the increase of Ra and Ta. However, it is found that the value of ħ tends to increase with the increase of Ra, but it decreases with the increase of Ta.Table 2
Part of the valid values of parametersℏ and γ used in the subsequent discussion.
Rayleigh number Ra
The valid values of parameters(ℏ,γ)
Ta
=
0
Ta
=
0.5
Ta
=
1.0
Ta
=
2.0
Ta
=
4.0
500
(
-
0.40
,
6
)
(
-
0.45
,
6
)
(
-
0.40
,
7
)
(
-
0.45
,
7
)
(
-
0.50
,
7
)
1000
(
-
0.20
,
6
)
(
-
0.30
,
8
)
(
-
0.35
,
9
)
(
-
0.35
,
9
)
(
-
0.40
,
9
)
2000
(
-
0.20
,
8
)
(
-
0.25
,
9
)
(
-
0.3
,
10
)
(
-
0.35
,
11
)
(
-
0.40
,
11
)


### 4.2. The Effect of Rayleigh Number and Taylor Number on Heat Transfer Rate

Once the values ofħ and γ are determined, the Nusselt number, temperature distribution, and velocity field at any spatial point only depend on Ra and Ta. In this work, the values of Ta are limited to a interval [0,4], while the values of Ra are taken to be as high as 2000. The justification of these values can refer to previous studies [1, 19, 64] on rotating effect or high-Rayleigh-number effect on thermal convection in porous media.The Nusselt number Nu as a function of Taylor number Ta for three differentRa is shown in Figure 5. First, it is clearly seen that, for each curve, the Nusselt number decreases with the increase of Taylor number, indicating that the presence of rotation can reduce the heat transfer rate in porous media at all considered Ra. Second, it is shown in Figure 5 that, at given values of Ta, the reducing rate of Nu is nearly the same for the three curves. Finally, it can be observed that the reducing effect of rotation on the Nu is more obvious at smaller Ta. In fact, the Nu tends to have an asymptotic value for very large Ta.Figure 5
Nusselt number Nu as a function ofTa for three different Ra.The Nusselt number Nu as a function of Rayleigh numberRa for three different Ta is illustrated in Figure 6. It can be seen that the Nu increases rapidly with the increase of Ra. This increment is very reasonable because higher Ra indicates larger temperature difference which serves as the source of thermal buoyancy force. In contrast to the effect of Ta, there exists no asymptotic value of Nu for large Ra. Furthermore, the comparison of three curves shown in Figure 6 can also confirm the reducing effect of rotation on Nu.Figure 6
Nusselt number Nu as a function ofRa for three different Ta.


### 4.3. The Effect of Rayleigh Number and Taylor Number on Temperature Distributions and Convective Flows

The effects ofRa and Ta on temperature distributions are illustrated in Figures 7 and 8. It can be seen that the temperature decays with the z for both groups of parameters. Further, the smaller the Ta or the larger the Ra is, the more rapidly the temperature decreases. It is also obvious that the wall-normal temperature gradient ∂θ/∂z|z=0 is higher at the smaller Ta or the larger Ra, which is in close agreement with the results shown in Figures 5 and 6.Figure 7
The effect ofRa on temperature distribution at r=1 with Ta=2.Figure 8
The effect of Ta on temperature distribution atr=1 with Ra=1000.Figure9 shows the effects of Ta on the velocity components w. It can be seen that the velocities oscillate with z. Furthermore, for a larger Ta, the w has one smaller crest and one larger trough, while, for smaller Ta, the w has one larger crest and one smaller trough. This indicates that rotation tends to give rise to downward flow and reduce the amplitude of upward thermal convection. Hence, the thermal convection with stronger rotation is weaker, which is consistent with the effect of rotation on heat transfer of thermal convection.Figure 9
The effect of Ta on the velocity componentw at r=0.15 with Ra=1000.The effects of Rayleigh number on velocity profiles are displayed in Figure10. With the increase of Ra, the amplitude of upward velocity increases as expected. Besides, it can been seen from Figure 10 that, for larger Ra, the interval of z corresponding to the positive value of w is greater than that for smaller Ra, indicating that thermal convection with larger Ra penetrates farther away from the bottom.Figure 10
The effect ofRa on the velocity component w at r=0.15 with Ta=2.Finally, the pattern of stream line is illustrated in Figure11 where the directions of the arrows coincide with the directions of flows in velocity fields. It can be seen that the convection rolls arise near the bottom boundary. However, the upward thermal convection in the present study arises from the interaction between the destabilizing effect of thermal buoyancy force and the suppressing effect of rotation. Therefore, the increase of Ta would reduce the upward thermal convection in a rotating porous half space.Figure 11
Stream lines forRa=1400 and Ta=1.It is noted that this work mainly concentrates on the steady sate thermal convection in a porous media; that is, all the physical variables are independent of time. Therefore, after dimensionalizing the governing equations (1) the Darcy-Prandtl number (defined as the ratio of Prandtl number and Darcy number [59]) is not included in the dimensionless equations (3)–(5). Besides, the typical values of Darcy-Prandtl number in traditional porous media applications are quite large [22, 59], a fact which provides the justification for taking the infinite Darcy-Prandtl number limit in some studies [1, 16]. However, for the general case of unsteady convection, the Darcy-Prandtl number does not vanish in the governing equation. It is interesting that the homotopy analysis method can be used even for Darcy-Prandtl number. In this case, the time derivative operator related to Darcy-Prandtl number in the momentum should be chosen as the auxiliary linear operator of homotopy analysis. The initial guess solutions are the functions of time which satisfy the auxiliary linear operator, and the base functions should also be modified to include the time component. Then, using the procedure of homotopy analysis method, the approximate solutions for the case of Darcy-Prandtl number can be obtained.It is also expected that the homotopy analysis method can be applied for nonlinear stability analysis of thermal convection in porous media. The objective of the nonlinear analysis is to provide quantitative results regarding the time evolution of convection amplitude. The expansion of eigenfunctions for the nonlinear governing equation in terms of a small perturbation parameter [13, 22, 23] and the truncated Galerkin representation of eigenfunctions for the nonlinear governing equation with a minimal double Fourier series [24, 59, 65] are the most widely used methods to perform nonlinear stability analysis. The resulting equations describing the amplitude of convection for both methods are nonlinear autonomous system of ordinary differential equations with respect to time. Arriving at this step, the homotopy analysis method can be applied to solve the sets ordinary differential equations. By choosing proper base function and auxiliary linear operator, the homotopy analysis method can provide analytical solutions to the nonlinear equations in terms of base function.


## 4.1. The Discussion on the Parametersħ and γ

The discussed problem is governed by coupling nonlinear partial differential equations (3)–(5), which cannot be analytically solved by traditional methods. In this work, analytical solutions are presented by virtue of the homotopy analysis method (HAM) based on the homotopy, a basic concept of topology. Unlike previous analytic techniques, the HAM provides us with a family of solution series by introducing the nonzero auxiliary parameter ħ. The parameter ħ plays an important role in the HAM. However, it cannot be chosen arbitrary because the convergence of solution series is dependent on its value.In order to clarify how to determine the valid value or interval of parameterħ, the so-called ħ-curves in the HAM are used, in which some related physical quantities versus parameter ħ are plotted. According to Liao [50, 62] who first proposed the HAM, if there exists a nearly horizontal line segment in the ħ-curve, the horizontal line segment corresponds to the valid interval of ħ, and the solution series are convergent for each ħ in the valid interval.In the present study, the curve Nu versusħ serves as the ħ-curve. Our strategy for determining the valid value of ħ proceeds as follows. Firstly, for a set of given parameters Ra,Ta, and γ, we plot the Nu-ħ curve. Secondly, we check whether there is a nearly horizontal line segment in the Nu-ħ curve. If there is none, the choice of γ is not appropriate; we change the value of γ and replot the Nu-ħ curve until the horizontal line segment appears. Finally, the valid interval of ħ and the value of γ are determined accordingly. Repeating this procedure for different Ra and Ta, we can determine all the valid values of ħ and γ used in our computation.Using the procedure mentioned above, the Nu-ħ curve with Ra=1000 and Ta=4 is shown in Figure 2. It can be seen that, for γ=9, there exist nearly horizontal line segments, and the valid interval of ħ is about [-0.3,-0.5]. For each value of ħ in the interval [-0.3,-0.5], the solution series (20) and (27) of the problem are expected to be convergent. However, it is found that, for other values of γ, no such horizontal line segment exists, indicating that, in those cases, the solution series are divergent for any ħ. Figure 3 shows the effect of valid ħ on the temperature distribution, where three values of ħ in the valid interval are taken. It can be seen that the difference among three curves is very small, which implies that if ħ is chosen from the valid interval its effect on the solution series is negligible.Figure 2
The display of Nu-ħ curve which is used to determine the valid interval of ħ.Figure 3
The effect of validħ on the temperature distribution.It is noteworthy that, in practical applications, the infinite solution series are truncated when the subscript satisfiesm>N (N is a positive integer). The accuracy of the truncation scheme is checked by choosing sufficiently large N such that the solutions do not change significantly when N is replaced by N+2. The Nu-ħ curves with Ra=1000, Ta=4, and γ=9 for different truncated order N are shown in Figure 4. It can be seen that the 11th-order truncated approximation has already been satisfactory.Figure 4
The Nu-ħ curves for different truncated order N with Ra=1000, Ta=4, and γ=9.In Table1, the values of Nusselt number with Ra=1000 and Ta=4 are calculated for different valid parameter ħ and truncated order N. From this table, it can be seen that the convergence rate of solution series for different valid ħ is different. Therefore, the parameter ħ can not only control whether the solution series are convergent but also adjust the convergence rate of solution series. Furthermore, it is easily obtained that the relative error between 11th-order truncated approximation and 9th-order truncated approximation is smaller than 1, which is accurate enough to predict the heat transfer rate.Table 1
The values of Nusselt number Nu for different truncated orderN and valid parameter ℏ of solution series with Ra=1000 and Ta=4.
Truncated orderN
The values of Nusselt number Nu
ℏ
=
-
0.30
ℏ
=
-
0.35
ℏ
=
-
0.40
ℏ
=
-
0.45
7
3.73733
3.38959
3.27674
3.25194
8
3.46868
3.29490
3.27273
3.27276
9
3.37163
3.27839
3.26864
3.26225
10
3.34877
3.26358
3.26591
3.26846
11
3.27582
3.26235
3.26534
3.26538The valid values of parametersħ and γ used in the following computation are partly listed in Table 2. It can be seen that the value of γ tends to increase with the increase of Ra and Ta. However, it is found that the value of ħ tends to increase with the increase of Ra, but it decreases with the increase of Ta.Table 2
Part of the valid values of parametersℏ and γ used in the subsequent discussion.
Rayleigh number Ra
The valid values of parameters(ℏ,γ)
Ta
=
0
Ta
=
0.5
Ta
=
1.0
Ta
=
2.0
Ta
=
4.0
500
(
-
0.40
,
6
)
(
-
0.45
,
6
)
(
-
0.40
,
7
)
(
-
0.45
,
7
)
(
-
0.50
,
7
)
1000
(
-
0.20
,
6
)
(
-
0.30
,
8
)
(
-
0.35
,
9
)
(
-
0.35
,
9
)
(
-
0.40
,
9
)
2000
(
-
0.20
,
8
)
(
-
0.25
,
9
)
(
-
0.3
,
10
)
(
-
0.35
,
11
)
(
-
0.40
,
11
)


## 4.2. The Effect of Rayleigh Number and Taylor Number on Heat Transfer Rate

Once the values ofħ and γ are determined, the Nusselt number, temperature distribution, and velocity field at any spatial point only depend on Ra and Ta. In this work, the values of Ta are limited to a interval [0,4], while the values of Ra are taken to be as high as 2000. The justification of these values can refer to previous studies [1, 19, 64] on rotating effect or high-Rayleigh-number effect on thermal convection in porous media.The Nusselt number Nu as a function of Taylor number Ta for three differentRa is shown in Figure 5. First, it is clearly seen that, for each curve, the Nusselt number decreases with the increase of Taylor number, indicating that the presence of rotation can reduce the heat transfer rate in porous media at all considered Ra. Second, it is shown in Figure 5 that, at given values of Ta, the reducing rate of Nu is nearly the same for the three curves. Finally, it can be observed that the reducing effect of rotation on the Nu is more obvious at smaller Ta. In fact, the Nu tends to have an asymptotic value for very large Ta.Figure 5
Nusselt number Nu as a function ofTa for three different Ra.The Nusselt number Nu as a function of Rayleigh numberRa for three different Ta is illustrated in Figure 6. It can be seen that the Nu increases rapidly with the increase of Ra. This increment is very reasonable because higher Ra indicates larger temperature difference which serves as the source of thermal buoyancy force. In contrast to the effect of Ta, there exists no asymptotic value of Nu for large Ra. Furthermore, the comparison of three curves shown in Figure 6 can also confirm the reducing effect of rotation on Nu.Figure 6
Nusselt number Nu as a function ofRa for three different Ta.


## 4.3. The Effect of Rayleigh Number and Taylor Number on Temperature Distributions and Convective Flows

The effects ofRa and Ta on temperature distributions are illustrated in Figures 7 and 8. It can be seen that the temperature decays with the z for both groups of parameters. Further, the smaller the Ta or the larger the Ra is, the more rapidly the temperature decreases. It is also obvious that the wall-normal temperature gradient ∂θ/∂z|z=0 is higher at the smaller Ta or the larger Ra, which is in close agreement with the results shown in Figures 5 and 6.Figure 7
The effect ofRa on temperature distribution at r=1 with Ta=2.Figure 8
The effect of Ta on temperature distribution atr=1 with Ra=1000.Figure9 shows the effects of Ta on the velocity components w. It can be seen that the velocities oscillate with z. Furthermore, for a larger Ta, the w has one smaller crest and one larger trough, while, for smaller Ta, the w has one larger crest and one smaller trough. This indicates that rotation tends to give rise to downward flow and reduce the amplitude of upward thermal convection. Hence, the thermal convection with stronger rotation is weaker, which is consistent with the effect of rotation on heat transfer of thermal convection.Figure 9
The effect of Ta on the velocity componentw at r=0.15 with Ra=1000.The effects of Rayleigh number on velocity profiles are displayed in Figure10. With the increase of Ra, the amplitude of upward velocity increases as expected. Besides, it can been seen from Figure 10 that, for larger Ra, the interval of z corresponding to the positive value of w is greater than that for smaller Ra, indicating that thermal convection with larger Ra penetrates farther away from the bottom.Figure 10
The effect ofRa on the velocity component w at r=0.15 with Ta=2.Finally, the pattern of stream line is illustrated in Figure11 where the directions of the arrows coincide with the directions of flows in velocity fields. It can be seen that the convection rolls arise near the bottom boundary. However, the upward thermal convection in the present study arises from the interaction between the destabilizing effect of thermal buoyancy force and the suppressing effect of rotation. Therefore, the increase of Ta would reduce the upward thermal convection in a rotating porous half space.Figure 11
Stream lines forRa=1400 and Ta=1.It is noted that this work mainly concentrates on the steady sate thermal convection in a porous media; that is, all the physical variables are independent of time. Therefore, after dimensionalizing the governing equations (1) the Darcy-Prandtl number (defined as the ratio of Prandtl number and Darcy number [59]) is not included in the dimensionless equations (3)–(5). Besides, the typical values of Darcy-Prandtl number in traditional porous media applications are quite large [22, 59], a fact which provides the justification for taking the infinite Darcy-Prandtl number limit in some studies [1, 16]. However, for the general case of unsteady convection, the Darcy-Prandtl number does not vanish in the governing equation. It is interesting that the homotopy analysis method can be used even for Darcy-Prandtl number. In this case, the time derivative operator related to Darcy-Prandtl number in the momentum should be chosen as the auxiliary linear operator of homotopy analysis. The initial guess solutions are the functions of time which satisfy the auxiliary linear operator, and the base functions should also be modified to include the time component. Then, using the procedure of homotopy analysis method, the approximate solutions for the case of Darcy-Prandtl number can be obtained.It is also expected that the homotopy analysis method can be applied for nonlinear stability analysis of thermal convection in porous media. The objective of the nonlinear analysis is to provide quantitative results regarding the time evolution of convection amplitude. The expansion of eigenfunctions for the nonlinear governing equation in terms of a small perturbation parameter [13, 22, 23] and the truncated Galerkin representation of eigenfunctions for the nonlinear governing equation with a minimal double Fourier series [24, 59, 65] are the most widely used methods to perform nonlinear stability analysis. The resulting equations describing the amplitude of convection for both methods are nonlinear autonomous system of ordinary differential equations with respect to time. Arriving at this step, the homotopy analysis method can be applied to solve the sets ordinary differential equations. By choosing proper base function and auxiliary linear operator, the homotopy analysis method can provide analytical solutions to the nonlinear equations in terms of base function.


## 5. Conclusion

This paper deals with steady state thermal convection in a rotating porous half space. By using the homotopy analysis method (HAM), the coupled nonlinear momentum and energy equations are solved to obtain the analytical solutions for velocity and temperature fields of thermal convection. The convergence and validity of HAM is guaranteed by plotting the so-calledħ-curves. The effects of Taylor number and Rayleigh number on the heat transfer, convective flow, and temperature distribution are investigated by using the truncated solution series. The main results show that as the Rayleigh number increases the heat transfer rate is enhanced rapidly. However, with the increase of Taylor number, the heat transfer rate is reduced, and the Nusselt number tends to have an asymptotic value for very large Taylor number.The obtained result in this paper can convey some physical meanings of thermal convection in the considered system. The convective flow when a fluid is heated from below is manifest as thermal convection. Since the gravity buoyancy in terms of Rayleigh number is the fundamental force to give rise to thermal convection, the convective flow with higher Rayleigh number (larger temperature difference) is more intense than that with relatively lower Rayleigh number. Furthermore, the temperature in the vertical directions decreases more rapidly with the increase of Rayleigh number, indicating the increasing heat transfer rate. Finally, for a fluid-saturated half space, the rotation can give rise to downward flow which is opposite to upward flow of thermal convection. Consequently, when the fluids are subject to an external constraint like rotation, thermal convection in porous media can be suppressed and the heat transfer rate can be reduced accordingly.



---

*Source: 905458-2015-03-09.xml*