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

Authors: Rajesh Sharma; Anuar Ishak
Journal: Mathematical Problems in Engineering
 (2013)

---

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


---

## Body

## 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 Figure1. The fluid motion arises due to the motion of the plate. The x-coordinate is measured along the plate from its leading edge, and y-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 x-direction; hence, all physical quantities will be independent of x. 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
(1)∂v∂y=0.Momentum Conservation. One has(2)∂u∂t+v∂u∂y=ν∂2u∂y2-(νεkp+Fε2kp1/2u)u+βg(h-h∞).Energy Conservation. One has
(3)∂h∂t+v∂h∂y=α∂2h∂y2+1cp(νεkpu2+Fε2kp1/2u3+ν(∂u∂y)2).
The corresponding boundary conditions on the vertical surface and in the free stream can be defined now as follows:
(4)t<0:u=0,h=h∞∀y,t≥0:Aty=0:u=up,v=-v0,h=hw,y⟶∞:u=0,h⟶h∞,
where u and v are the velocity components along the x and y axes, h is the temperature inside the boundary layer, t is time and ρ, ν, kp, ε, F, β, cp, and g 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), up, v0, hw, and h∞ 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)U=uV0,V=vV0,Y=V0yν,T=V02tν,Up=upV0,θ=h-h∞hw-h∞.
The resulting equations are
(6)∂U∂T-∂U∂Y=∂2U∂Y2-(1Da+FrDaU)U+Grθ∂θ∂T-∂θ∂Y=1Pr∂2θ∂Y2+Ec[1DaU2+FrDaU3+(∂U∂Y)2],
where Ec=V02/cp(hw-h∞) is the Eckert number, Pr=ν/α is the Prandtl number, Da=kpε/L2 is the Darcy number, Gr=νgβ(hw-h∞)/V03 is Grashof number, Fr=Fεkp1/2/L is Forchheimer inertial porous media number, T is non-dimensional time, U is dimensionless velocity, and θ is dimensionless temperature.The corresponding transformed temporal and spatial boundary conditions are(7)T<0:U=0,θ=0∀Y,T≥0:AtY=0:U=Up,θ=1,Y⟶∞:U=U∞⟶0,θ⟶0.
Using Fourier’s law, the heat flux at the plate surface may be written as follows:
(8)qw=-k0∂h∂y|y=0=-k0υ0ρμ(hw-h∞)θ′(0),
where k0 is the coefficient of thermal conductivity.The heat transfer coefficient is given by(9)hf=qw(hw-h∞)=-k0υ0ρμθ′(0).
The local Nusselt number/heat transfer rate can be written as
(10)Nux=hfxk0=-Rexθ′(0).


## 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 U¯ and θ¯, which are assumed to be known values of the functions U 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 (Ye,Ye+1) is given by
(11)∫YeYe+1w1{∂U∂T-∂U∂Y-∂2U∂Y2+(1DaRex2+FrDaRexU)U-Grθ}dY=0,∫YeYe+1w2{∂θ∂T-∂θ∂Y-1Pr∂2θ∂Y2-Ec[1DaRex2U2+FrDaRexU3+(∂U∂Y)2]}dY=0,
where w1 and w2 are arbitrary test functions and may be viewed as the variation in U and θ, respectively. After reducing the order of integration and nonlinearity, we arrive at the following system of equations:
(12)∫0Y∞{w1∂U∂T-w1∂U∂Y+∂w1∂Y∂U∂Y+w1(1DaRex2+FrDaRexU-)U-Grw1θ}dY-(w1∂U∂Y)Ye+1Ye=0,∫YeYe+1{{[1DaRex2U-U+FrDaRexU-2U+∂U-∂Y∂U∂Y]}w2∂θ∂T-w2∂θ∂Y+1Pr∂w2∂u∂θ∂Y-Ecw2[1DaRex2U-U+FrDaRexU-2U+∂U-∂Y∂U∂Y]}dY-(w2Pr∂θ∂Y)Ye+1Ye=0.


### 4.2. Finite Element Formulation

The finite element model may be obtained from (11) by substituting finite element approximations of the form
(13)U=∑j=12Ujeψje,θ=∑j=12θjeψje
with w1=w2=ψie(i=1,2), where Uje and θje are the velocity and temperature, respectively, at jth node of typical eth element (Ye,Ye+1) and ψie are the shape functions for this typical element (Ye,Ye+1) and are taken as
(14)ψ1e=Ye+1-YYe+1-Ye,ψ2e=Y-YeYe+1-YeYe≤Y≤Ye+1.
The finite element model of the equations for the element thus formed is given by
(15)[[K11][K12][K21][K22]][[U][θ]]+[[M11][M12][M21][M22]][[U˙][θ˙]]=[[F1][F2]],
where [Klm], [Mlm], and [Fl] (l,m=1,2) are the matrices of order 2×2, 2×2 and 2×1, respectively. [U] and [θ] are matrices of order 2×1. Here, U˙ and θ˙ are derivatives of Uand θ with respect to T. All these matrices may be defined as follows:
(16)(K11)IJ=∫YeYe+1(-ψIe∂ψJe∂Y+∂ψIe∂Y∂ψJe∂Y+ψIe(1DaRex2+FrDaRexU-)ψJe)dY(K12)IJ=-Gr∫YeYe+1(ψIeψJe)dY,(K21)IJ=-Ec∫YeYe+1ψIe(1DaRex2U-ψJe+FrDaRexU-2ψJe+∂U-∂Y∂ψJe∂Y)dY,(K22)IJ=∫0Y∞(-ψIe∂ψJe∂Y+∂ψIe∂Y1Pr∂ψJe∂Y)dY,(M11)IJ=∫0Y∞(ψIeψJe)dY,(M12)IJ=(M21)IJ=0,(M22)IJ=∫0Y∞(ψIeψJe)dY,(F1)I=(ψIe∂U∂Y)YeYe+1,(F2)I=1PrψIe∂θ∂Y,
where
(17)U¯=∑I=12UI¯ψIe.It has been found that, for large value ofY∞(≥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, Y∞ 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.


## 4.1. Variational Formulation

The variational formulation associated with (6) over a typical two node linear elements (Ye,Ye+1) is given by
(11)∫YeYe+1w1{∂U∂T-∂U∂Y-∂2U∂Y2+(1DaRex2+FrDaRexU)U-Grθ}dY=0,∫YeYe+1w2{∂θ∂T-∂θ∂Y-1Pr∂2θ∂Y2-Ec[1DaRex2U2+FrDaRexU3+(∂U∂Y)2]}dY=0,
where w1 and w2 are arbitrary test functions and may be viewed as the variation in U and θ, respectively. After reducing the order of integration and nonlinearity, we arrive at the following system of equations:
(12)∫0Y∞{w1∂U∂T-w1∂U∂Y+∂w1∂Y∂U∂Y+w1(1DaRex2+FrDaRexU-)U-Grw1θ}dY-(w1∂U∂Y)Ye+1Ye=0,∫YeYe+1{{[1DaRex2U-U+FrDaRexU-2U+∂U-∂Y∂U∂Y]}w2∂θ∂T-w2∂θ∂Y+1Pr∂w2∂u∂θ∂Y-Ecw2[1DaRex2U-U+FrDaRexU-2U+∂U-∂Y∂U∂Y]}dY-(w2Pr∂θ∂Y)Ye+1Ye=0.


## 4.2. Finite Element Formulation

The finite element model may be obtained from (11) by substituting finite element approximations of the form
(13)U=∑j=12Ujeψje,θ=∑j=12θjeψje
with w1=w2=ψie(i=1,2), where Uje and θje are the velocity and temperature, respectively, at jth node of typical eth element (Ye,Ye+1) and ψie are the shape functions for this typical element (Ye,Ye+1) and are taken as
(14)ψ1e=Ye+1-YYe+1-Ye,ψ2e=Y-YeYe+1-YeYe≤Y≤Ye+1.
The finite element model of the equations for the element thus formed is given by
(15)[[K11][K12][K21][K22]][[U][θ]]+[[M11][M12][M21][M22]][[U˙][θ˙]]=[[F1][F2]],
where [Klm], [Mlm], and [Fl] (l,m=1,2) are the matrices of order 2×2, 2×2 and 2×1, respectively. [U] and [θ] are matrices of order 2×1. Here, U˙ and θ˙ are derivatives of Uand θ with respect to T. All these matrices may be defined as follows:
(16)(K11)IJ=∫YeYe+1(-ψIe∂ψJe∂Y+∂ψIe∂Y∂ψJe∂Y+ψIe(1DaRex2+FrDaRexU-)ψJe)dY(K12)IJ=-Gr∫YeYe+1(ψIeψJe)dY,(K21)IJ=-Ec∫YeYe+1ψIe(1DaRex2U-ψJe+FrDaRexU-2ψJe+∂U-∂Y∂ψJe∂Y)dY,(K22)IJ=∫0Y∞(-ψIe∂ψJe∂Y+∂ψIe∂Y1Pr∂ψJe∂Y)dY,(M11)IJ=∫0Y∞(ψIeψJe)dY,(M12)IJ=(M21)IJ=0,(M22)IJ=∫0Y∞(ψIeψJe)dY,(F1)I=(ψIe∂U∂Y)YeYe+1,(F2)I=1PrψIe∂θ∂Y,
where
(17)U¯=∑I=12UI¯ψIe.It has been found that, for large value ofY∞(≥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, Y∞ 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 Figures2–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, Rex=1.0, and Up=0.5.Figure 2
3D visualization of nondimensional velocity profile at different time(T).Figure 3
3D visualization of nondimensional temperature profile at different time(T).Figure 4
Velocity distribution at different time(T).Figure 5
Temperature distribution at different time(T).Figure 6
Velocity distribution for differentUp.Figure 7
Temperature distribution for differentUp.Figure 8
Velocity distribution for differentGr.Figure 9
Temperature distribution for differentGr.Figure 10
Velocity distribution for differentDa.Figure 11
Temperature distribution for differentDa.Figure 12
Velocity distribution for differentFr.Figure 13
Temperature distribution for differentFr.Figure 14
Temperature distribution for differentEc.Figure 15
Nusselt number for differentUp.Figure 16
Nusselt number for differentGr.Figure 17
Nusselt number for differentDa.Figure 18
Nusselt number for differentPr.Figure 19
Nusselt number for differentEc.The three-dimensional visualization of velocity(U) and temperature (θ) profiles with both space (Y-coordinate) and time (T) is shown in Figures 2 and 3, where time T varies from 1 to 6. The corresponding two-dimensional plots of velocity U and temperature θ with space coordinate Y for different value of time T is shown in Figures 4 and 5. It is observed that there is a small change in U and θ after T=5 for all computations. Thus, the results for T=5 are essentially the steady-state values. The velocity U starts from 0.5 at Y=0 and attains a maximum near the plate, thereafter decay back to zero asymptotically at Y=8.The influence of plate velocity(Up) on the velocity and temperature with Ycoordinate 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 Y 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 Y→∞.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 Y 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 Y 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 alongY 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 (Y=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 (Up) 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 Table1. From the simulation, it can be observed that the results are in good agreement for the dimensionless velocity (U) and dimensionless temperature (θ).Table 1
Comparison of FEM and FDM computation.
Y
U
(
Y
)
θ
(
Y
)
FDM
FEM
FDM
FEM
0.1
0.5734
0.5737
0.9282
0.9286
1.1
0.3750
0.3753
0.5588
0.5591
2.1
0.2111
0.2113
0.3123
0.3125
3.1
0.1090
0.1091
0.1618
0.1619
4.1
0.0523
0.0524
0.0781
0.0781
5.1
0.0231
0.0231
0.0347
0.0347
6.1
0.0090
0.0090
0.0136
0.0136
7.1
0.0024
0.0024
0.0038
0.0038


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



---

*Source: 371971-2013-09-19.xml*