Numerical Simulation of Transient Free Convection Flow and Heat Transfer in a Porous Medium

Abstract

The coupled momentum and heat transfer in unsteady, incompressible flow along a semi-infinite vertical porous moving plate adjacent to an isotropic porous medium with viscous dissipation effect are investigated. The Darcy-Forchheimer nonlinear drag force model which includes the effects of inertia drag forces is employed. The governing differential equations of the problem are transformed into a system of nondimensional differential equations which are solved numerically by the finite element method (FEM). The non-dimensional velocity and temperature profiles are presented for the influence of Darcy number, Forchheimer number, Grashof number, Eckert number, Prandtl number, plate velocity, and time. The Nusselt number is also evaluated and compared with finite difference method (FDM), which shows excellent agreement.

1. Introduction

Transport modeling in porous media finds numerous applications in diverse branches of physics, engineering, and applied mathematics. These include chemical reaction processes [1], heat transfer in petroleum reservoirs [2], reentry astronautical flows [3], geothermal energy simulations [4, 5], biothermomechanics of tissue [6], flame dynamics and ignition flows in porous media combustion systems [7, 8], thermal insulation systems [9], and solar energy collector systems [10]. Generally, the Darcy model is used to model the fluid flow through a porous medium, which is valid under conditions of low velocities and small porosity. This model is the most widely adopted in mathematical modeling of porous media hydrodynamics. Recently, Motsa and Shateyi [11] examined the problem of magnetomicropolar fluid flow, heat, and mass transfer with suction through a porous medium under the effect of chemical reaction, Hall, and ion slip currents. Non-Darcian flows however are important in high-velocity applications and have been considered by many researchers [1, 12–15] for a diverse range of multiphysical flows in porous media.

Unsteady transport in porous media has also attracted considerable attention. Various studies have employed these models. Bhargava et al. [16] investigated the oscillatory reactive hydromagnetic natural convection boundary layers in porous media with cross flow gradient effects. Zueco et al. [17] studied transient Darcy-Forchheimer flow in a rotating porous medium using network simulation. Toki [18] studied the transient natural convection flow near an oscillating porous infinite vertical plate (or wall) with heating of the plate using a Laplace transform technique. Raptis and Perdikis [19] investigated analytically the transient convection-radiation flow in Darcian regime. Bég et al. [20] used network simulation to model the hydromagnetic convection in a porous medium channel. Cheng and Lin [21] considered thermal dispersion effects on transient convection using the Darcy-Forchheimer model. El-Kabeir et al. [22] studied numerically the magnetic field effects on a translating sheet embedded in porous media. Graham and Higdon [23] used a finite element method to analyze oscillatory forcing of flow through porous media with the time-dependent Navier-Stokes equations. Cheng et al. [24] used a boundary element method (BEM) to model unsteady flow in 3-dimensional reservoir bounded porous medium with various boundary conditions.

The previously mentioned studies ignored the effect of viscous dissipation, which plays an important role in the dynamics of fluids. The heat generated by viscous dissipation produces an increase in temperature near the walls with a consequent decrease of the viscosity and a strong stratification in the viscosity profile, which affects the heat transfer rate. The modeling of viscous dissipation in a porous medium is completely different from that of pure fluid flow [25, 26]. However, there are some studies [27–29] that modeled it similar to that in pure fluid. In the literature, the study of viscous dissipation effect on fluid flow through a porous medium was done by Nield [30] and Al-Hadhrami et al. [31]. Nield [25] commented that the model proposed by Al-Hadhrami et al. [31] is probably adequate for most practical purposes.

In this paper, finite element method (FEM) is utilized to simulate the unsteady free convection flow of a viscous fluid past a semi-infinite, steadily moving porous plate adjacent to an isotropic, homogenous Non-Darcy porous medium with constant suction velocity normal to the plate. Viscous dissipation is included which is modeled in accordance with Al-Hadhrami et al. [31].

2. Mathematical Model

Consider an unsteady free convection flow of Newtonian, viscous fluid in a Darcy-Forchheimmer porous medium along a vertical semi-infinite moving permeable flat plate in the presence of viscous heating, which is sketched in Figure 1. The fluid motion arises due to the motion of the plate. The 𝑥-coordinate is measured along the plate from its leading edge, and 𝑦-coordinate is taken normal to it. In the analysis, the flow is laminar, unsteady, and two-dimensional; the porous medium is isotropic and homogeneous; the fluid and porous medium are everywhere in local thermodynamic equilibrium. The plate is considered semi-infinite in the 𝑥-direction; hence, all physical quantities will be independent of 𝑥. With these assumptions, the governing equations for free convective flow through Darcy-Forchheimmer porous medium describing the conservation of mass, momentum and energy can be written as follows.

Figure 1 

Flow model and coordinate system.

Mass Conservation. One has 𝜕𝑣𝜕𝑦=0. (1)

Momentum Conservation. One has𝜕𝑢𝜕𝑡+𝑣𝜕𝑢𝜕𝜕𝑦=𝜈²𝑢𝜕𝑦²−𝜈𝜀𝑘_𝑝+𝐹𝜀²𝑘_𝑝^1/2𝑢𝑢+𝛽𝑔ℎ−ℎ_∞. (2)

Energy Conservation. One has 𝜕ℎ𝜕𝑡+𝑣𝜕ℎ𝜕𝜕𝑦=𝛼²ℎ𝜕𝑦²+1𝑐_𝑝𝜈𝜀𝑘_𝑝𝑢²+𝐹𝜀²𝑘_𝑝^1/2𝑢³+𝜈𝜕𝑢𝜕𝑦². (3) The corresponding boundary conditions on the vertical surface and in the free stream can be defined now as follows: 𝑡<0:𝑢=0,ℎ=ℎ_∞∀𝑦,𝑡≥0:At𝑦=0:𝑢=𝑢_𝑝,𝑣=−𝑣₀,ℎ=ℎ_𝑤,𝑦⟶∞:𝑢=0,ℎ⟶ℎ_∞, (4) where 𝑢 and 𝑣 are the velocity components along the 𝑥 and 𝑦 axes, ℎ is the temperature inside the boundary layer, 𝑡 is time and 𝜌, 𝜈, 𝑘_𝑝, 𝜀, 𝐹, 𝛽, 𝑐_𝑝, and 𝑔 are the density, kinematic viscosity, permeability, porosity, the Forchheimer empirical constant, thermal expansion coefficient, specific heat at constant pressure, and the acceleration due to gravity, respectively. The middle three terms on the right-hand side of the energy equation (3) signify viscous dissipation effect in Darcy-Forchheimer media which are modeled according to Al-Hadhrami et al. [25]. In (4), 𝑢_𝑝, 𝑣₀, ℎ_𝑤, and ℎ_∞ are the plate velocity, suction, wall temperature, and free-stream temperature, respectively.

3. Transformation of Model

The partial differential equations (1)–(3) with boundary condition (4) are transformed into non-dimensional form by means of following dimensionless variables: 𝑢𝑈=𝑉₀𝑣,𝑉=𝑉₀𝑉,𝑌=⁰𝑦𝜈,𝑉𝑇=²⁰𝑡𝜈,𝑈_𝑝=𝑢_𝑝𝑉₀,𝜃=ℎ−ℎ_∞ℎ_𝑤−ℎ_∞. (5) The resulting equations are 𝜕𝑈−𝜕𝑇𝜕𝑈=𝜕𝜕𝑌²𝑈𝜕𝑌²−1+DaFr𝑈Da𝑈+Gr𝜃𝜕𝜃−𝜕𝑇𝜕𝜃=1𝜕𝑌𝜕Pr²𝜃𝜕𝑌²1+Ec𝑈Da²+Fr𝑈Da³+𝜕𝑈𝜕𝑌², (6) where Ec=𝑉²₀/𝑐_𝑝(ℎ_𝑤−ℎ_∞) is the Eckert number, Pr=𝜈/𝛼 is the Prandtl number, Da=𝑘_𝑝𝜀/𝐿² is the Darcy number, Gr=𝜈𝑔𝛽(ℎ_𝑤−ℎ_∞)/𝑉³₀ is Grashof number, Fr=𝐹𝜀𝑘_𝑝^1/2/L is Forchheimer inertial porous media number, 𝑇 is non-dimensional time, 𝑈 is dimensionless velocity, and 𝜃 is dimensionless temperature.

The corresponding transformed temporal and spatial boundary conditions are 𝑇<0:𝑈=0,𝜃=0∀𝑌,𝑇≥0:At𝑌=0:𝑈=𝑈_𝑝,𝜃=1,𝑌⟶∞:𝑈=𝑈_∞⟶0,𝜃⟶0. (7) Using Fourier’s law, the heat flux at the plate surface may be written as follows: 𝑞_𝑤=−𝑘₀𝜕ℎ𝜕𝑦_𝑦=0𝑘=−⁰𝜐₀𝜌𝜇ℎ_𝑤−ℎ_∞𝜃^′(0), (8) where 𝑘₀ is the coefficient of thermal conductivity.

The heat transfer coefficient is given by ℎ_𝑓=𝑞_𝑤ℎ_𝑤−ℎ_∞𝑘=−⁰𝜐₀𝜌𝜇𝜃^′(0). (9) The local Nusselt number/heat transfer rate can be written as Nu_𝑥=ℎ_𝑓𝑥𝑘₀=−Re_𝑥𝜃^′(0). (10)

4. Numerical Solution

The finite element method (FEM) is employed to solve the two-point boundary value problem defined by (6) under conditions (7). FEM is the most versatile of all modern numerical methods and is succinctly described by Bathe [32], Reddy [33], Bhargava et al. [16], and Sharma et al. [34]. The uniform nodal points are used to discretize the spatial coordinate. An implicit Crank-Nicolson scheme which is robust and unconditionally stable is used to discretize the time derivative. The time step has been taken as 0.0001. An iterative scheme has been used to solve the resulting system of nonlinear equations at each time level. The system of equations is first linearized by incorporating the functions 𝑈 and 𝜃, which are assumed to be known values of the functions 𝑈 and 𝜃. After applying the given boundary conditions, a resulting system of equations is solved using Gauss-elimination method. This gives new values of unknowns for the next iteration. This process continues until the absolute differences between two successive iterations become less than the accuracy of 0.0001. The following steps are involved in the analysis.

4.1. Variational Formulation

The variational formulation associated with (6) over a typical two node linear elements (𝑌_𝑒,𝑌_𝑒+1) is given by ^𝑌𝑒+1_𝑌𝑒𝑤₁𝜕𝑈−𝜕𝑇𝜕𝑈−𝜕𝜕𝑌²𝑈𝜕𝑌²+1DaRe²_𝑥+FrDaRe_𝑥𝑈𝑈−Gr𝜃𝑑𝑌=0,_{𝑌𝑒+1𝑌𝑒}𝑤₂𝜕𝜃−𝜕𝑇𝜕𝜃−1𝜕𝑌𝜕Pr²𝜃𝜕𝑌²1−EcDaRe²_𝑥𝑈²+FrDaRe_𝑥𝑈³+𝜕𝑈𝜕𝑌²𝑑𝑌=0, (11) where 𝑤₁ and 𝑤₂ are arbitrary test functions and may be viewed as the variation in 𝑈 and 𝜃, respectively. After reducing the order of integration and nonlinearity, we arrive at the following system of equations: ^𝑌∞₀𝑤₁𝜕𝑈𝜕𝑇−𝑤₁𝜕𝑈+𝜕𝑌𝜕𝑤₁𝜕𝑌𝜕𝑈𝜕𝑌+𝑤₁1DaRe²_𝑥+FrDaRe_𝑥𝑈𝑈−Gr𝑤₁𝜃−𝑤𝑑𝑌₁𝜕𝑈𝜕𝑌^𝑌𝑒_𝑌𝑒+1=0,_{𝑌𝑒+1𝑌𝑒}𝑤₂𝜕𝜃𝜕𝑇−𝑤₂𝜕𝜃+1𝜕𝑌Pr𝜕𝑤₂𝜕𝑢𝜕𝜃𝜕𝑌−Ec𝑤₂1DaRe²_𝑥𝑈𝑈+FrDaRe_𝑥𝑈²𝜕𝑈+𝑈𝜕𝑌𝜕𝑈−𝑤𝜕𝑌𝑑𝑌₂Pr𝜕𝜃𝜕𝑌^𝑌𝑒_𝑌𝑒+1=0. (12)

4.2. Finite Element Formulation

The finite element model may be obtained from (11) by substituting finite element approximations of the form 𝑈=²_𝑗=1𝑈^𝑒_𝑗𝜓^𝑒_𝑗,𝜃=²_𝑗=1𝜃^𝑒_𝑗𝜓^𝑒_𝑗 (13) with 𝑤₁=𝑤₂=𝜓^𝑒_𝑖(𝑖=1,2), where 𝑈^𝑒_𝑗 and 𝜃^𝑒_𝑗 are the velocity and temperature, respectively, at 𝑗th node of typical 𝑒th element (𝑌_𝑒,𝑌_𝑒+1) and 𝜓^𝑒_𝑖 are the shape functions for this typical element (𝑌_𝑒,𝑌_𝑒+1) and are taken as 𝜓^𝑒₁=𝑌_𝑒+1−𝑌𝑌_𝑒+1−𝑌_𝑒,𝜓^𝑒₂=𝑌−𝑌_𝑒𝑌_𝑒+1−𝑌_𝑒𝑌_𝑒≤𝑌≤𝑌_𝑒+1. (14) The finite element model of the equations for the element thus formed is given by 𝐾₁₁𝐾₁₂𝐾₂₁𝐾₂₂[𝑈][𝜃]+𝑀₁₁𝑀₁₂𝑀₂₁𝑀₂₂𝑈̇𝜃̇=𝐹₁𝐹₂, (15) where [𝐾_𝑙𝑚], [𝑀_𝑙𝑚], and [𝐹_𝑙] (𝑙,𝑚=1,2) are the matrices of order 2×2, 2×2 and 2×1, respectively. [𝑈] and [𝜃] are matrices of order 2×1. Here, 𝑈̇ and 𝜃̇ are derivatives of 𝑈and 𝜃 with respect to 𝑇. All these matrices may be defined as follows: 𝐾_11𝐼𝐽=^𝑌𝑒+1_𝑌𝑒−𝜓^𝑒_𝐼𝜕𝜓^𝑒_𝐽+𝜕𝑌𝜕𝜓^𝑒_𝐼𝜕𝑌𝜕𝜓^𝑒_𝐽𝜕𝑌+𝜓^𝑒_𝐼1DaRe²_𝑥+FrDaRe_𝑥𝑈𝜓^𝑒_𝐽𝐾𝑑𝑌_12𝐼𝐽=−Gr^𝑌𝑒+1_𝑌𝑒𝜓^𝑒_𝐼𝜓^𝑒_𝐽𝐾𝑑𝑌,_21𝐼𝐽=−Ec^𝑌𝑒+1_𝑌𝑒𝜓^𝑒_𝐼1DaRe²_𝑥𝑈𝜓^𝑒_𝐽+FrDaRe_𝑥𝑈²𝜓^𝑒_𝐽+𝜕𝑈𝜕𝑌𝜕𝜓^𝑒_𝐽𝐾𝜕𝑌𝑑𝑌,_22𝐼𝐽=^𝑌∞₀−𝜓^𝑒_𝐼𝜕𝜓^𝑒_𝐽+𝜕𝑌𝜕𝜓^𝑒_𝐼1𝜕𝑌Pr𝜕𝜓^𝑒_𝐽𝑀𝜕𝑌𝑑𝑌,_11𝐼𝐽=^𝑌∞₀𝜓^𝑒_𝐼𝜓^𝑒_𝐽𝑀𝑑𝑌,_12𝐼𝐽=𝑀_21𝐼𝐽𝑀=0,_22𝐼𝐽=^𝑌∞₀𝜓^𝑒_𝐼𝜓^𝑒_𝐽𝐹𝑑𝑌,_1𝐼=𝜓^𝑒_𝐼𝜕𝑈𝜕𝑌^𝑌𝑒+1_𝑌𝑒,𝐹_2𝐼=1𝜓Pr^𝑒𝐼𝜕𝜃,𝜕𝑌 (16) where 𝑈=²_𝐼=1𝑈_𝐼𝜓^𝑒_𝐼. (17)

It has been found that, for large value of 𝑌_∞(≥8), there is no appreciable change in the results up to desired order of accuracy; therefore, for the computational purpose and without loss of generality, 𝑌_∞ has been fixed as 8. The adequacy of the elements number is verified by comparing the results of different set of elements; 160 elements are found to be satisfactory. Therefore, the domain is represented by a set of 160 uniformly elements.

5. Result and Discussion

Computations have been carried out for various values of parameters such as plate velocity, Grashof number, Eckert number, Prandtl number, Darcy number, and Forchheimer number. The numerical results so obtained are plotted in Figures 2–19 for velocity, temperature, and Nusselt number. The default values of parameters are taken as Da=1.0, Gr=5.0, Fr=5.0, Ec=0.005, Pr=0.71, Re_𝑥=1.0, and 𝑈_𝑝=0.5.

Figure 2 

3D visualization of nondimensional velocity profile at different time (𝑇).

Figure 3 

3D visualization of nondimensional temperature profile at different time (𝑇).

Figure 4 

Velocity distribution at different time (𝑇).

Figure 5 

Temperature distribution at different time (𝑇).

Figure 6 

Velocity distribution for different 𝑈𝑝.

Figure 7 

Temperature distribution for different 𝑈𝑝.

Figure 8 

Velocity distribution for different Gr.

Figure 9 

Temperature distribution for different Gr.

Figure 10 

Velocity distribution for different Da.

Figure 11 

Temperature distribution for different Da.

Figure 12 

Velocity distribution for different Fr.

Figure 13 

Temperature distribution for different Fr.

Figure 14 

Temperature distribution for different Ec.

Figure 15 

Nusselt number for different 𝑈𝑝.

Figure 16 

Nusselt number for different Gr.

Figure 17 

Nusselt number for different Da.

Figure 18 

Nusselt number for different Pr.

Figure 19 

Nusselt number for different Ec.

The three-dimensional visualization of velocity (𝑈) and temperature (𝜃) profiles with both space (𝑌-coordinate) and time (𝑇) is shown in Figures 2 and 3, where time 𝑇 varies from 1 to 6. The corresponding two-dimensional plots of velocity 𝑈 and temperature 𝜃 with space coordinate 𝑌 for different value of time 𝑇 is shown in Figures 4 and 5. It is observed that there is a small change in 𝑈 and 𝜃 after 𝑇=5 for all computations. Thus, the results for 𝑇=5 are essentially the steady-state values. The velocity 𝑈 starts from 0.5 at 𝑌=0 and attains a maximum near the plate, thereafter decay back to zero asymptotically at 𝑌=8.

The influence of plate velocity (𝑈_𝑝) on the velocity and temperature with 𝑌coordinate is depicted in Figures 6 and 7. An increase in plate velocity consequently increases the fluid velocity as well as velocity gradient near the plate wall. Because of this steep velocity gradient, internal energy of fluid at the wall is enhanced in terms of kinetic energy and hence the temperature increases as shown in Figure 7.

The influence of the convection parameter or Grashof number (Gr) on velocity and temperature with 𝑌 direction is depicted in Figures 8 and 9. As expected, it is observed that an increase in Gr leads to a rise in the velocity due to enhancement in buoyancy force. In addition, graph shows that, for Gr≤3.0, velocity continuously decreases, while for Gr>3.0, the velocity increases near the wall and then reduces to free stream velocity. Figure 9 shows the influence of Grashof number on temperature function 𝜃. Higher Gr value means that more heat is removed, so that temperature must decrease which is clear from the graphs also. All profiles decay asymptotically to zero as 𝑌→∞.

The influence of Darcy number (Da) on the velocity profile for the case of non-Darcian convection with viscous heating effect is shown in Figure 10. As Da increases, a significant increase in velocity occurs. It shows that with greater permeability, there will be a reduction in the Darcian drag acting on the flow regime. This will serve to progressively accelerate the flow, which is consistent with the trend obtained in Figure 10. It is also observed that, for Da<1.0, velocity continuously decreases throughout the boundary layer region, and for Da≥1.0, it first increases near the plate and then decreases to zero in the free stream.

The variation of temperature (𝜃) with 𝑌 for various Darcy number (Da) is illustrated in Figure 11. From Figure 11, it is seen that temperature decreases moderately as Da increases. Increase in permeability reduces the presence of solid particles in the medium, thereby reducing the conduction heat transfer which in turn reduces the temperature.

The effect of the Forchheimer parameter (inertial drag parameter) Fr on velocity along 𝑌 direction is depicted in Figure 12. Since Fr represents the inertial drag, an increase in the Forchheimer parameter increases the resistance to the flow in the porous medium, which will serve to decelerate the flow, that is, reduce velocities throughout the boundary layer. Here, Fr=0 represents the case, when the flow is Darcian; that is, inertial effects are neglected. Thus, the velocity is found maximum in this case.

The variation of the temperature along 𝑌 for various values of Fr is shown in Figure 13. A rise in Forchheimer number corresponds to a rise in temperature. Maximum temperature therefore is associated with the maximum value of Fr. Since the fluid is decelerated with the increase in Fr, the fluid energy is dissipated as heat and serves to enhance temperature in the porous medium, throughout the domain. Hence, the results are consistent with the physics behind the concept of Forchheimer parameter. All profiles decay continuously from the maximum at the surface (𝑌=0) to zero far away from the surface.

The effect of Eckert number (Ec) on the temperature is illustrated in Figure 14. Eckert number represents the viscous dissipation effect. It is observed that temperature increases with the increase in viscous dissipation effect. Further, it can be seen that temperature is higher with viscous dissipation as compared to the case of no viscous dissipation, that is, Ec=0. The increase in the fluid temperature due to viscous dissipation is observed to be more pronounced for a small change in Ec. For Ec≤0.005, temperature continuously decreases in the boundary layer region, and, for value greater than 0.005, it first increases near the plate and then starts decreasing.

The rate of heat transfer/Nusselt number is governed by the value of −𝜃^′(0) and the results are shown in Figures 15, 16, 17, 18, and 19. It is observed that heat transfer increases with the increase in convection parameter (Gr), Prandtl number (Pr), and Darcy parameter (Da). The increase in the values of Nusselt number shows that the heat is transferred from the plate to the fluid; that is, a cooling of the plate is taking place. Thus, it can be seen that the parameters Gr, Pr, and Da can be effectively used in controlling temperature in many industrial problems. A decrease in the heat transfer is seen with the increase in plate velocity (𝑈_𝑝) and Eckert number (Ec).

In order to check the accuracy of the present numerical computations, the problem is also solved using the finite difference method and a comparison of the results is presented in Table 1. From the simulation, it can be observed that the results are in good agreement for the dimensionless velocity (𝑈) and dimensionless temperature (𝜃).

Table 1 

Comparison of FEM and FDM computation.

6. Conclusions

A mathematical model has been presented for transient natural convection flow from a translating plate adjacent to a non-Darcian porous medium with viscous heating effects. The nondimensionalized conservation equations have been solved under physically realistic boundary conditions using the finite element method (FEM). On the basis of present simulations, the following conclusion can be made. (i)The velocity increases with the increase of plate velocity, Darcy number, and Grashof number, but it decreases with the increase of Forchheimer number. (ii)The fluid temperature increases with the increase of plate velocity and Eckert number, as a result of this Nusselt number decreases; that is, rate of heat transfer decreases. Therefore, the viscous heating effect may be used to reduce the rate of cooling in, for example, materials processing and nuclear heat transfer systems. (iii)The convection parameter, Prandtl number, can be used effectively for increasing the rate of heat transfer as required in many applications like nuclear reactor, aerodynamics, and so forth, where it is required to control the enormous temperature. (iv)Increase in Eckert number (Ec) heats the porous regime; that is, increases temperature, which result in a decrease in Nusselt number. Thus, it may be used to reduce the rate of cooling. But for Ec>0.005, fluid temperature near the wall is predicted to exceed the wall temperature inferring that the direction of heat transfer is reversed from fluid to the sheet.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors wish to express their thanks to the anonymous referees for their valuable comments and suggestions. The financial support received from the Universiti Kebangsaan Malaysia (Project code DIP-2012-31) is gratefully acknowledged.