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

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.


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][13][14][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][28][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].

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.

Mass Conservation. One has
V  = 0. (1) Momentum Conservation.One has Energy Conservation.One has The corresponding boundary conditions on the vertical surface and in the free stream can be defined now as follows: where  and V 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),   , V 0 , ℎ  , and ℎ ∞ are the plate velocity, suction, wall temperature, and free-stream temperature, respectively.

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: The resulting equations are where /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 Using Fourier's law, the heat flux at the plate surface may be written as follows: where  0 is the coefficient of thermal conductivity.
The heat transfer coefficient is given by The local Nusselt number/heat transfer rate can be written as

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

Variational Formulation.
The variational formulation associated with (6) over a typical two node linear elements (  ,  +1 ) is given by where  1 and  2 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: 4.2.Finite Element Formulation.The finite element model may be obtained from (11) by substituting finite element approximations of the form with , 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 The finite element model of the equations for the element thus formed is given by where where 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.

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.
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 ().

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.

Figure 1 :
Figure 1: Flow model and coordinate system.

Table 1 :
Comparison of FEM and FDM computation.