Numerical simulation of thermal conve*ction of viscoelastic fluids in an open-top porous medium with constant heat flux

NIU Jun （牛駿）, SHI Zai-hong （石在虹）, TAN Wen-chang （譚文長）

1. Petroleum Exploration and Production Research Institute of SINOPEC, Beijing 100083, China, E-mail: niujun.syky@sinopec.com

2. State Key Laboratory for Turbulence and Complex Systems and Department of Mechanics and Aerospace Engineering, College of Engineering, Peking University, Beijing 100871, China

Numerical simulation of thermal conve*ction of viscoelastic fluids in an open-top porous medium with constant heat flux

NIU Jun （牛駿）1, SHI Zai-hong （石在虹）1, TAN Wen-chang （譚文長）2

1. Petroleum Exploration and Production Research Institute of SINOPEC, Beijing 100083, China, E-mail: niujun.syky@sinopec.com

2. State Key Laboratory for Turbulence and Complex Systems and Department of Mechanics and Aerospace Engineering, College of Engineering, Peking University, Beijing 100871, China

(Received July 23, 2013, Revised April 8, 2014)

This paper makes a numerical study of the buoyancy-driven convection of a viscoelastic fluid saturated in an open-top porous square box under the constant heat flux boundary condition. The effects of the relaxation and retardation times on the onset of the oscillatory convection, the convection heat transfer rate and the flow pattern transition are investigated. It is shown that a large relaxation time can destabilize the fluid flow leading to an early onset of the thermal convection and a high heat transfer rate, while a large retardation time tends to stabilize the flow and suppress the convection onset and the heat transfer. After the convection sets in, the flow bifurcation appears earlier with the increase of the relaxation time and the decrease of the retardation time, resulting in more complicated flow patterns in the porous medium. Furthermore, with the increase of the ratio of the relaxation time to the retardation time, the fluid may be blocked from flowing through the open-top boundary, which may be caused by the viscoelastic effect. Finally, the comparison of our results with those under isothermal heating boundary conditions reveals that the heat transfer rate corresponding to a constant heat flux boundary is always higher.

thermal convection, viscoelastic fluids, porous medium, constant heat flux

Introduction

It is of a great importance to study the thermal convection of viscoelastic fluids in porous media, as its applications are found in many engineering fields, like the paper and textile coating, the composite manufacturing processes, and the bio-engineering[1-4]. In the oil reservoir engineering, the heavy oil is also found to exhibit viscoelastic rheological characteristics[5]. However, the thermal convection of viscoelastic fluids in porous media is far less fruitfully studied than that of the Newtonian fluids[6-8]. The main reason may be due to the lack of a simple model for the description of the viscoelastic flow behavior in porous media as the Darcy’s law for the Newtonian fluids.

Recently, the modified Darcy’s law[9-12]has drawn much attention in the research of thermal instability problems of viscoelastic fluid saturated porous media heated from below. This modified Darcy’s law is a macroscopic phenomenological model, which incorporates the viscoelastic effects suggested by Alisaev and Mirzadjanzade[13]and is based on the rheological behavior of viscoelastic liquids. With the modified Darcy’s law, Niu et al.[14]conducted stability analyses of viscoelastic fluids in square porous media under three different heating boundary conditions, which shows that an oscillatory convection mode always sets in before a steady one. Kim et al.[9]conducted linear and non-linear stability analyses of the thermal convection in an Oldroyd-B fluid saturated porous layer with isothermal bottom heating. However, these studies are based on a linear or weakly non-linear analysis, where the Darcy-Rayleigh number should be around the critical values for the onset of the convection.

When the Rayleigh number gets higher and the analytical results become invalid, numerical simulations must be conducted in order to investigate theflow and heat transfer process. Fu et al.[15]numerically analyzed the thermal convection of a viscoelastic fluid in an impermeable porous square box heated from below with a constant bottom temperature. It is found that after the steady convection sets in, the oscillatory convection tends to be suppressed and the interaction between the two types of thermal convections makes the flow patterns in the porous layer very complicated. Furthermore, the thermal convection of viscoelastic fuids in porous layers would undergo earlier bifurcations with the increase of the Darcy-Rayleigh number than that of the Newtonian fuids and the occurrence of the bifurcation is earlier for a larger ratio of the relaxation time to the retardation time. As a consequence, the heat transfer characteristics are quite different from those of the Newtonian fluids. The heat transfer rate (as represented by the Nusselt number) is essentially higher for a larger ratio of the relaxation time to the retardation time, which obeys different scaling laws delimited by the viscoelastic parameters and the Darcy-Rayleigh number. The thermal convection for a viscoelastic fluid saturated in an open-top (permeable) porous medium with an isothermal heating boundary[16]was also studied, which shows that the onset of the thermal convection occurs earlier than that under a closed-top condition, as the open-top boundary imposes a weaker constraint against the fluid flow inside the porous medium. Therefore, the heat transfer rate after the onset of the thermal convection is larger compared to the closed-top case. However, to the best of our knowledge, no results were obtained for the thermal convection of a viscoelastic fluid in a porous medium under a constant heat flux boundary condition.

To study the thermal convection in viscoelastic fluid saturated open-top porous media under a constant heat flux boundary condition, a numerical simulation was made in our previous letter[17], where the viscoelastic effects on the convection evolution were briefly discussed.

The objective of this paper is to systematically analyze the entire thermal convection process, beginning from the convection onset to the convection bifurcation at high Rayleigh numbers. Therefore, a linear stability analysis is first conducted to study the onset of the thermal convection in a viscoelastic fluid saturated open-top porous media under the constant heat flux boundary condition. Then, a numerical algorithm is built and the corresponding calculation program is developed for analyzing the convection heat transfer and the bifurcation. The influence of different heating conditions is studied by comparing our results with those in the previous work.

1. Model formation

We consider a model of a bounded two-dimensional rectangular porous medium of length a and height H*, with a permeable isothermal top boundary at temperature0T*, where the pressure at the top boundary is assumed as constant[16]. The impermeable bottom boundary is under a constant heat flux q. The two vertical boundaries are adiabatic and impermeable. The porous medium has a permeability K, which is saturated by an incompressible viscoelastic fluid with a constant dynamic viscosity μ, a coefficient of thermal expansion β and a density ρ. The fluid-saturated porous medium has a thermal conductivity k and a thermal diffusivity κ. From the modified Darcy’s law, the governing equations of this model are[14-16]

where v*=(u*,w*)is the Darcy velocity, p*the pressure, g the gravitational acceleration,ε and λ are, respectively, the strain retardation time and the stress relaxation time, z is the unit vector along the z- direction which is vertically upward, and ρ0the density at temperature0T. The boundary conditions for this model are

The constant top pressure boundary condition (5b) is used to describe the open top and is converted into an equivalent velocity boundary condition (1+ε?/?t*)· u*(x*, H*)=0 by substituting (5b) into (2).

Our aim is to analyze the whole heat transfer process in the viscoelastic fluid saturated porous mediumwith an open top and a constant bottom heat flux. This process can be divided into two parts: the onset of the thermal convection and the convection evolution after the onset[9,14,15]. A linear stability analysis is carried out in Section 3 for solving the onset of the thermal convection, while a numerical algorithm based on the finite difference method is adopted in Section 3 for calculating the convection evolution.

2. Linear stability analysis

In order to perform the linear stability analysis, Eqs.(1)-(4) are first non-dimensionalized by scaling length with H*, time with2/Hκ*, temperature with qH*/k and velocity with κ/H*. Then θ and w are introduced to represent, respectively, the infinitesimal disturbances of the dimensionless temperature and the vertical velocity over the pure conduction solution. As a result, Eqs.(1)-(4) are simplified as

The Rayleigh number for this model is Ra=Kβρ0g·2/ H qkμ κ. The non-dimensionalized viscoelastic parameters are ε=εκ/H2and λ=λκ/H2. The boundary conditions are

where =/aaH**is the dimensinless length of the porous medium.

In the normal mode analysis, the temperature and velocity disturbances are supposed to be horizontally periodic. Under the adiabatic vertical boundary conditions, θ and w take the forms

where σ is the temporal growth rate and α=mπ/ a , m is an arbitrary integer. ()zΘ is a function of z only, and =/Dz?? is a differential operator.

Substituting Eqs.(8), (9) into (6) yields

Equation (11) has the solution in the form

where Aj(j=1,2,3,4) are unknown constants, kjare the roots of the following equation

According to the boundary conditions (8a)-(8d),2=A 1A,43=AA, and1A,3A,1k,3k satisfy the relations

As1A,3A cannot both be zero, the linear Eqs.(14) and (15) have non-zero solutions, which requires the coefficient determinant to be zero

The onset of the thermal convection is a steady convection as that for the Newtonian fluid if =0σ. But for the viscoelastic fluid, the oscillatory convection, which is more of interest, may appear earlier. For investigating the onset of the oscillatory convection, we set σ as a pure imaginary number iω and solve Eq.(16) numerically.

3. Numerical algorithm

The linear stability analysis performed in Section 3 can only be used to study the thermal convection near the onset point. After that, the non-linear effect may not be negligible, so the governing equationshave to be solved numerically for further investigations. However, due to the complexity of the governing equations, no simulation software is available for calculating the thermal convection of viscoelastic fluids in a porous medium.

Here, we build a numerical algorithm based on the finite difference method for simulating the thermal convection evolution.

Introducing a stream function ψ to represent the Darcy velocity, the dimensionless governing equations are

The boundary conditions of the current problem are given as

The bulk-averaged Nusselt number Nu is used to estimate the convection heat transport characteristics in the porous medium, which is given by

The angle brackets indicate the long-time average.

The finite difference algorithm is adopted in the current study. The energy Eq.(18) is solved using the alternating direction implicit (ADI) method as

where

and n denotes the nth time step.

The updated temperature distribution after a time step tΔ is obtained from Eq.(22), then Eq.(17) becomes a first-order ordinary differential equation with respect to ξ=▽2ψ, which can be solved with a second-order Runge-Kutta method[18]

The stream function ψ is calculated with a fast Fourier transform (FFT) based method[19], which is used later to obtain the velocities u and w. The solution domain is discretized into NN× uniform subdomains, where N is equal to a certain power of 2 in order for the FFT-based method to be applicable.

Based on the above algorithm, a calculation program is developed for the convection simulation.

4. Results and discussions

For investigating the effects of the parameters λ and ε on the onset of the thermal convection and the convection evolution characteristics of a viscoelastic fluid in an open-top porous layer under the bottom constant heat flux, three set values are selected as Case 1: =0.3λ, =0.2ε, Case 2: =0.2λ, =0.1ε,Case 3: =0.3λ, =0.1ε. As Cases 2 and 3 have the same retardation time value, the results can be used to analyze the effect of the relaxation time. Similarly, the comparison of the results of Cases 1 and 3 may reveal the effect of the retardation time.

4.1 Onset of oscillatory convection

The onset of the oscillatory convection is marked by a critical Rayleigh number cRa which is calculated from Eq.(16) over all modes m with fixed a, λ and ε. ThecRa for Cases 1-3 as a function of a along with the preferred mode m is shown in Fig.1. It is indicated that with the increase of a, the critical Rayleigh number decreases with the periodic vibrations. The global minimum Rayleigh numbers for the three cases are 25.63, 24.32 and 17.58, respectively.

Fig.1 The critical Rayleigh numbers for the onset of oscillatory convection for Cases 1-3 as functions of the porous medium length

Figures 1(a) and 1(c) show the variation ofcRa for Cases 1 and 3 with the same λ. For any fixed a thecRa in Fig.1(a) is always larger than that in Fig.1(c), which means that the increase of the retardation time ε can stabilize the flow and postpone the onset of the oscillatory convection. Comparing the variation ofcRa in Figs.1(b) and 1(c) for Cases 1 and 3 with the same ε, it can be seen that the increase of the relaxation time can disturb the flow, leading to an earlier onset of the thermal convection.

4.2 Convection heat transfer

After the onset of the oscillatory convection, our numerical algorithm is used to solve Eqs.(17)-(20) for Ra up to 400 to study the convection evolution. There, the thermal convection in a square porous medium with =1a is calculated as a representative case. In order to ensure the accuracy of our calculations, we choose =128N for 100Ra≤ and =256N for Ra＞100. Furthermore, the calculations are for Ra in the increasing direction and the increase intervals of Ra are in the range from 2 to 25. The calculation begins from the critical Rayleigh number for the three cases determined in Section 5.1. The input data at each Ra are taken from the results of the previous Ra.

Fig.2 Comparison of the heat transfer curves for the three cases of the viscoelastic fluid, showing the effect of the viscoasticity on the heat transfer rate. Also shown is the curve for a Newtonian fluid

The calculated results for the bulk-averaged Nusselt number Nu varying with the Darcy-Rayleigh number Ra up to 400 are shown in Fig.2 for the three set values of λ and ε as well as the Newtonian fluid case. For Cases 1 and 2, interestingly, parts of the two heat transfer curves almost overlap that for the Newtonian fluid. It can be inferred from this phenomenon that the steady convection after onset can suppress the oscillatory convection and dominate the convection process, as observed in the previous work[15]. This trend persists for Ra up to about 45. With the further increase of Ra, the heat transfer curve for Case 2 separates first from the curve for the Newtonian fluid and is located above the latter. The heat transfer curve for Case 1 will also separate from that for the Newtonian fluid when Ra is larger than about 50. However, unlike the curve for Case 2, it is located below the curve for the Newtonian fluid when Ra is in the range between 50 and 175. There is a sharp increase in the heat transfer rate for Case 1 at Ra=185, such that this curve jumps up and intersects the curve for the Newtonian fluid at Ra around 200. The heat transfer curve for Case 3, however, has nopart overlapping that for the Newtonian fluid and is located totally above other three curves shown in Fig.2. These distinct properties of the heat transfer process for Cases 1, 2 and 3 are caused by a series of transitions of the flow pattern inside the porous medium, which in turn are affected by the viscoelasticity of the fluid, as will be discussed in detail in the subsequent section.

Comparing the curves in Fig.1 for Cases 1 and 3, in which =0.3λ and ε is equal to 0.2, 0.1, it can be seen that a larger value of ε corresponds to a smaller Nu. Therefore, the increase of the retardation time tends to suppress the heat transfer. The curves for Cases 2 and 3, in which =0.1ε and λ is equal to 0.2, 0.3, show that the increase of the relaxation time can enhance the heat transfer rate.

Fig.3 The time history of Nu and contours of the stream function for Case 1 at =45Ra

Fig.4 Successive snapshots of the stream function contours for Case 1 at =125Ra with a time interval of 0.2, showing that the flow is in a one-cell roll alternating with a twocell roll pattern

Fig.5 The time history of Nu and power spectrum for Case 1 at =125Ra

4.3 Flow pattern transition and bifurcation

For the viscoelastic fluid, due to the effect of viscoelasticity, the oscillatory convection mode and the steady convection mode may co-exist after the onset of the thermal convection[17]. For Case 1, the steady convection mode dominates the heat transfer process for Ra up to 50. In this range, the flow exhibits a steady one-cell roll pattern and the heat transfer curve overlaps that for the Newtonian fluid. Figure 3(a) shows the time history of Nu for Case 1 at =45Ra. Because the input data at each Ra value are taken from the results for the previous Ra value in our calculation, Nu varies with time but, as can be seen from the figure that it eventually reaches a constant at large enough t. The constant Nu is a sign of the steady convection and is equal to the value for the Newtonian fluid at =45Ra. The corresponding flow pattern can be indicated by the stream function ψ, which is plotted in Fig.3(b). A one-cell roll flow pattern is clearly seen and the outer contours of ψ vertically cross the open-top boundary.

When Ra is greater than 50, the heat transfer curve for Case 1 deviates from and becomes lower than that for the Newtonian fluid. The convection can no longer remain in a steady one-cell roll state, but ex-hibits a flow pattern of a one-cell roll alternating with a two-cell roll. Figures 4(a)-4(d) show the snapshots of the stream function contours for Case 1 at =Ra 125, with a time interval equal to 0.2. A one-cell roll alternating with a two-cell roll flow pattern is clearly seen and the outer contours of the stream function do not cross the open-top boundary vertically, but bends towards the center of the flow circulation. From Eq.(8), the horizontal component of the flow velocity u at the top boundary is zero if ε is equal to zero. That is why the thermal convection for the Newtonian fluid flows vertically across the top boundary if it is permeable. For the viscoelastic fluid, however, u at the top boundary is a function of timeand ε, and thus can be nonzero due to the viscoelasticity. Figure 5(a) shows the time history of Nu for Case 1 at Ra=125, which reveals that the thermal convection is in a quasi-periodic state. Plotted in Fig.5(b) is the corresponding power spectrum of Nu, in which three prevailing frequencies can be identified: f1=1.33, f2=3.1 and f3=f1+f2=4.43, f1= 1.33. On the other hand, the flow pattern for the Newtonian fluid at =75Ra is a steady two-cell roll (it is a one-cell roll for =50Ra), resulting in a quick increase of the heat transfer rate. Therefore, though the oscillatory convection due to the viscoelasticity of the fluid affects significantly the heat transfer process, the heat transfer rate for Case 1 can be lower than that for the Newtonian fluid at the same value of Ra.

Fig.6 The time history of Nu and contours of the stream function for Case 1 at Ra=185

The heat transfer curve for Case 1 jumps sharply up to almost merge with that for the Newtonian fluid at Ra=185, which indicates the occurrence of another flow pattern transition. Shown in Fig.6(a) is the time history of Nu for Case 1 at Ra=185. In essence, this curve approaches a constant value after a long enough time and its characteristics are different from those shown in Fig.5(a) though tiny oscillations can be found inherent in it. Fig.6(b) shows a snapshot of the contours of the stream function for Case 1 at =Ra 185. A two-cell roll flow pattern is clearly seen and, unlike that shown in Fig.4, this flow pattern persists at any instant and represents a typical flow pattern. Furthermore, it is interesting to find that the contours of the stream function now change to cross the top boundary vertically, similar to the case for a Newtonian fluid in a steady two-cell roll state[20]. Therefore, it may be concluded that the steady convection re-dominates the thermal convection process after the transition, but the effect of the viscoelasticity is not completely suppressed, which is manifested by the tiny oscillations in the Nu vs. t curve and which may explain why the bulk-averaged heat transfer rate is slightly higher than that for the Newtonian fluid.

The two-cell roll flow pattern for Case 1 persists for Ra up to about 275, where the thermal convection experiences the third bifurcation and the flow pattern changes to an unsteady two-cell roll. This unsteady two-cell roll pattern further enhances the heat transfer rate, which increases the distance between the heat transfer curves for Case 1 and that of the Newtonian case as Ra increases. Figure 7(a) shows the time history of Nu for Case 1 at =325Ra, in which the Nusselt number vibrates quasi-periodically with time. The corresponding power spectrum is plotted in Fig.7(b), which reveals that there are also three prevailing frequencies in the wave envelope:1=1.54f, f2=2.01 and f3=3.81. These prevailing frequencies are different from those in Fig.5. It seems that the heat transfer rate depends on the flow pattern inside the porous medium, which in turn is affected by the viscoelasticity of the fluid.

Fig.7 The time history of Nu and power spectrum for Case 1 at Ra=325

The thermal convection for Case 2 is also dominated by the steady convection after the onset for Ra up to approximately 45. Therefore, the flow exhibits a steady one-cell roll pattern and the corresponding heat transfer curve in this range follows that for the Newtonian fluid. When Ra is further increased, however, the flow pattern for Case 2 will transit into an unsteady one-cell roll alternating with a two-cell roll, and the heat transfer curve will deviate from and become higher than that for the Newtonian fluid. This situation lasts for Ra up to 400. The time histories of Nu for Case 2 at Ra equal to 150, 250 and 350 are shown in Figs.8(a), 8(c) and 8(e), respectively. Similar quasi-periodic waveform of the Nu is found except that the period of the waveform in Fig. 8(a) is a bit larger than that in Figs.8(c) and 8(e). The corresponding power spectra of the Nu vs. t curves are shown in Figs.8(b), 8(d) and 89(f), respectively, and it is shown that the prevailing frequency1f for Fig. 8(a) is equal to 0.5 while it is equal to 0.6 for both Fig.8(c) andFig.8(e).

Fig.8 The time histories of Nu ((a), (c) and (e)) for Case 2 at Ra equal to 150, 250 and 350, respectively. The corresponding power spectrum shown ((b), (d) and (f)) reveals only one prevailing frequency in this case

Fig.9 Snapshots of the stream function contours for Case 3 at Ra=50, 100 and 250, showing the flow patterns after each bifurcation

The thermal convection for Case 3 exhibits very different phenomena from those for Cases 1 and 2. A thermal convection bifurcation occurs for Ra between 40 and 45, resulting in a transition of the flow pattern from the one-cell roll to an unsteady two-cell roll. Another bifurcation occurs for Ra between 65 and 75 and the flow pattern transits from the unsteady two-cell roll to an unsteady three-cell roll after the bifurcation. The third bifurcation appears for Ra between 200 and 225, and the flow pattern varies from the unsteady three-cell roll to an unsteady three-cell roll alternating with a four-cell roll. Figures 9(a)-9(c) show the snapshots of the stream function contours for Case 3 at =Ra50, 100 and 250. The flow patterns seem to be under a close-top boundary similar to that found in the previous work[15], which may be induced by the high elastic effect.

From above analyses, it can be concluded that the increase of the relaxation time and the decrease of the retardation time can facilitate earlier thermal convection bifurcation and complicate the flow pattern.

Fig.10 Comparison of the heat transfer rates of the viscoelastic fluids with isothermal heating boundary and constant heat flux boundary

4.4 Comparison with viscoelastic fluid in an open-top porous medium with isothermal bottom heating

To see the influence of different bottom heating conditions, the heat transfer curves are compared in Fig.10 with those with the isothermal heating boundary. It can be seen that for each combination of λ and ε the heat transfer rate corresponding to the constant heat flux boundary is higher than that corresponding to the isothermal heating boundary at any fixed Ra. Our numerical results show that the onset of the thermal convection under the constant heat flux boundary condition is always earlier, which is consistent with what is observed in the previous work[14]. This earlier appearance of the convection will definitely enhance the heat transfer rate comparing to the conduction. Moreover, it is shown in Section 4.2 that each flow bi-furcation corresponding to the constant heat flux appears in a lower Rayleigh number than that with the isothermal bottom heating[16]. This earlier bifurcation tends to complicate the flow pattern, which may also lead to the enhancement of the heat transfer.

5. Conclusions

The thermal convection of a viscoelastic fluid saturated in an open-top porous square box with a constant heat flux is analyzed in this paper. The linear stability analysis is first performed for determining the critical Rayleigh number for the onset of the oscillatory convection. Then a finite difference algorithm is applied to numerically calculate the fluid flow and the heat transfer corresponding to different combinations of viscoelastic parameters and Rayleigh numbers. The results for three set values of the relaxation and retardation time constants of the viscoelastic fluid reveal some unique characteristics of the thermal convection process. First, the increase of the relaxation time can disturb the flow and result in an earlier onset of the thermal convection, while the increase of the retardation time would stabilize the flow and postpone the convection onset. Second, after the onset of the thermal convection the viscoelastic effect may either enhance or diminish the heat transfer rate comparing to the case for a Newtonian fluid at the same Rayleigh number, depending on the values of λ and ε. A larger λ tends to destabilize the thermal convection, leading to a higher heat transfer rate and earlier flow bifurcation, while a larger ε would stabilize the flow and simplify the flow pattern. Third, at the top boundary, the viscoelasticity of the fluid can result in a case with no flow across the boundary though it is permeable, as shown in Case 3. Finally, the comparison shows that the heat transfer rate of the viscoelastic fluid in the open-top porous layer under the constant heat flux boundary condition is always higher than that under the isothermal heating boundary condition at fixed Rayleigh number and viscoelastic parameters.

Our study may find applications in the oil thermal recovery engineering, such as the heavy oil recovery and the in-situ oil shale recovery, in which the porous reservoirs need to be heated. Our simulation algorithm can help predicting the temperature as well as the pressure distribution, and guide the site heating implementation.

This work focuses on the numerical simulation, and a corresponding experimental validation would yield further insight into the thermal convection in the viscoelastic fluid saturated porous media. But before that one needs to solve the problem we are facing regarding how to accurately measure the temperature distribution inside the porous media, which is also a cutting-edge topic.

[1] NIELD D. A., BEJAN A. Convection in porous media[M]. 4th Edition, New York, USA: Springer-Verlag, 2013.

[2] KHALED A. R. A., VAFAI K. The role of porous media in modeling flow and heat transfer in biological tissues[J]. International Journal of Heat and Mass Transfer, 2003, 46(26): 4989-5003.

[3] SHARMA Y. D., KUMAR V. Overstability analysis of thermo-bioconvection saturating a porous medium in a suspension of gyrotactic microorganisms[J]. Transport in Porous Media, 2011, 90(2): 673-685.

[4] HAYAT T., AWAIS M. and OBAIDAT S. Similar solution for three-dimensional flow in an Oldroyd-B fluid over a stretching surface[J]. International Journal for Numerical Methods in Fluids, 2012, 70(7): 851-859.

[5] MAKARYNSKA D., GUREVICH B. and BEHURA J. et al. Fluid substitution in rocks saturated with viscoelastic fluids[J]. Geophysics, 2010, 75(2): E115-E122.

[6] SHAN Lian-tao, TONG Dend-ke and XUE Li-li. Unsteady flow of non-Newtonian visco-elastic fluid in dual-porosity media with the fractional derivative[J]. Journal of Hydrodynamics, 2009, 21(5): 705-713.

[7] ZHANG Li-juan, YUE Xiang-an. Mechanism for viscoelastic polymer solution percolating through porous media[J]. Journal of Hydrodynamics, Ser. B, 2007, 19(2): 241-248.

[8] YANG Shu-ren, WANG Chun-sheng and CUI Hai-qing. Numerical silmulation of steady flow for viscoelastic fluid in an eccentric annulus with inner rod moving axially[J]. Journal of Hydrodynamics, Ser. B, 2005, 17(4): 514-518.

[9] KIM M. C., LEE S. B. and KIM S. J. et al. Thermal instability of viscoelastic fluids in porous media[J]. International Journal of Heat and Mass Transfer, 2003 46(26): 5065-5072.

[10] KHUZHAYOROV B., AURIAULT J. L. and ROYER P. Derivation of macroscopic filtration law for transient linear viscoelastic fluid flow in porous media[J]. International Journal of Engneering Science, 2000, 38(5): 487-504.

[11] BERTOLA B., CAFARO E. Thermal instability of viscoealstic fluids in horizontal porous layers as initial value problems[J]. International Journal of Heat and Mass Transfer, 2006, 49(21): 4003-4012.

[12] YOON D. Y., KIM M. C. and CHOI C. K. The onset of oscillatory convection in a horizontal porous layer saturated with viscoelastic liquid[J]. Transport in Porous Media, 2004, 55(3): 275-284.

[13] ALISHAEV M. G., MIRZADJANZADE A. K. For the calculation of delay phenomenon in filtration theory[J]. Izvestya Vuzov Neft i Gaz, 1975 6(1): 71-73.

[14] NIU J., FU C. and TAN W. Stability of thermal convection of an Oldroyd-B fluid in a porous medium with Newtonian heating[J]. Physics Letters A, 2010, 374(45): 4607-4613.

[15] FU C., ZHANG Z. and TAN W. Numerical simulation of thermal convection of an Oldroyd-B fluid in a porous square box heated from below[J]. Physics of Fluids, 2007, 19(10): 104107.

[16] NIU J., FU C. and TAN W. Thermal convection of a viscoelastic fluid in an open-top porous layer heated from below[J]. Journal of Non-Newtonian Fluid Mechanics, 2010, 165(5): 203-211.

[17] NIU Jun, SHI Zai-hong and TAN Wen-chang. The vis-coelastic effects on thermal convection of an Oldroyd-B fluid in open-top porous media[J]. Journal of Hydrodynamics, 2013, 25(4): 639-642.

[18] FU De-xun, MA Yan-wen. Computational fluid dynamics[M]. Beijing, China: Higher Education Press, 2002(in Chinese).

[19] MORGAN P. An introduction to fast fourier transform methods for partial differential equations, with applications[M]. Letchworth, UK: Research Studies Press, 1986.

[20] CHERKAUI A. S. M., WILCOCK W. S. D. Characteristics of high Rayleigh number two-dimensional convection in an open-top porous layer heated from below[J]. Journal of Fluid Mechanics, 1999, 394: 241-260.

10.1016/S1001-6058(15)60455-3

* Project supported by the National Key Basic Research Development Program of China (973 Program, Grant Nos. 2006CB705803, 2013CB531200).

Biography: NIU Jun (1985-), Male, Ph. D., Senior Engineer