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

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 analysismethod.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 theNusselt number decreases rapidlywith 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.


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][2][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 fluidsaturated 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 higherorder 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 Section 2, 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.

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  and porosity  is heated from below with constant temperature  0 and cooled from above at infinity with constant temperature  ∞ .Meanwhile, the porous medium rotates around -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  ≪ /Ω 2 ) from the axis of rotation, one can justify that the gravity buoyancy is dominant and centrifugal buoyancy can be neglected [22].
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] in which V, , and  are the Darcian velocity, pressure, and temperature. 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, where  0 is the density at temperature  0 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][6][7][8][9][10][11][12][13][14][15][16][17][18][19][21][22][23][24][25][26][27][28].
Choose /, /, and  0 − ∞ to scale the velocity, pressure, and temperature respectively; we can obtain from (1)-( 2) the dimensionless governing equations in the following form: where we have used symbol  = (− ∞ )/( 0 − ∞ ) 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 Previous studies have indicated that two-dimensional flow structures of thermal convection are possible in a fluid-saturated porous layer [16][17][18][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 (, ) = ((1/)(/), −(1/)(/)), where  and  are the components of V in the and -directions, respectively.Then, (3)-( 5) are reduced to the following form: The boundary conditions for ( 7) can be formulated as (i)  = 1 and  = 0, at  = 0; (ii)  = 0 and / = 0, at  → ∞.
By introducing the stream function, the considered problem is greatly simplified so that it can be solved analytically by the homotopy analysis method (HAM).

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:  = (, ) and  = (, ) where  = /.Then (7) are readily transformed into the following equations: The corresponding boundary conditions are (i)  = 1 and  = 0, at  = 0; (ii)  = 0 and / = 0, at  → ∞.

The Zero-Order Deformation Equations.
Following the basic idea of HAM, we assume that (, ) and (, ) can be expressed by a set of base functions {     − |  ≥ 0,  ≥ 0,  ≥ 0}: where   , and   , are coefficients and  is a spatial-scale parameter.
It can be observed that as  increases from 0 to 1, Φ(, ; ) and Θ(, ; ) vary continuously from the initial guess solutions  0 (, ) and  0 (, ) to the exact solutions (, ) and (, ) for the considered problem.This kind of continuous variation is called deformation in topology.Now, we expand Φ(, ; ) and Θ(, ; ) in a power series of the embedding parameter  as follows: where 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)

The High-Order Deformation Equations.
In (20), the functions   and   ( = 1, 2, 3, . ..) remain unknown.Differentiating the zero-order deformation equations ( 13)  times with respect to  and setting  = 0, the -order deformation equations can be obtained as subject to boundaries conditions where Let  *  and  *  denote particular solutions of (21); then we can express the general solutions as follows: where the coefficients  1 ,  2 ,  3 , and  4 can be determined by the boundary conditions (22).By virtue of symbolic computation software, one can compute   and   from ( 21)-( 24) order by order.

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 where ⟨⋅⟩  denotes area-averaged integration over a circular area A, | =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  and temperature difference Δ 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 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.It is noteworthy that, in practical applications, the infinite solution series are truncated when the subscript satisfies  >  ( is a positive integer).The accuracy of the truncation scheme is checked by choosing sufficiently large  such that the solutions do not change significantly when  is replaced by  + 2. The Nu-ℏ curves with Ra = 1000, Ta = 4, and  = 9 for different truncated order  are shown in Figure 4.It can be seen that the 11th-order truncated approximation has already been satisfactory.
In Table 1, the values of Nusselt number with Ra = 1000 and Ta = 4 are calculated for different valid parameter ℏ and truncated order .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.The 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.

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 different Ra 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.
The Nusselt number Nu as a function of Rayleigh number Ra 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.

The Effect of Rayleigh Number and Taylor Number on
Temperature Distributions and Convective Flows.The effects of Ra and Ta on temperature distributions are illustrated in Figures 7 and 8.It can be seen that the temperature decays with the  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 /| =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 9 shows the effects of Ta on the velocity components .It can be seen that the velocities oscillate with .Furthermore, for a larger Ta, the  has one smaller crest and one larger trough, while, for smaller Ta, the  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.
The effects of Rayleigh number on velocity profiles are displayed in Figure 10.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  corresponding to the positive value of  is greater than that for smaller Ra, indicating that thermal convection with larger Ra penetrates farther away from the bottom.Finally, the pattern of stream line is illustrated in Figure 11 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.
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.

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 fluidsaturated 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.

Figure 1 :
Figure 1: Schematic diagram of the physical problem.

9 Figure 2 : 4 Figure 3 :
Figure 2: The display of Nu-ℏ curve which is used to determine the valid interval of ℏ.

Figure 5 :
Figure 5: Nusselt number Nu as a function of Ta for three different Ra.

Figure 6 :
Figure 6: Nusselt number Nu as a function of Ra for three different Ta.

Figure 7 :
Figure 7: The effect of Ra on temperature distribution at  = 1 with Ta = 2.

Figure 10 :Figure 11 :
Figure 10: The effect of Ra on the velocity component  at  = 0.15 with Ta = 2.

Table 1 :
The values of Nusselt number Nu for different truncated order  and valid parameter ℎ of solution series with Ra = 1000 and Ta = 4.

Table 2 :
Part of the valid values of parameters ℎ and  used in the subsequent discussion.