The Mixed Convection Flow in a Fluid-saturated Non-Darcy Porous Medium through a Horizontal Channel

In this study the aim is to study the unsteady mixed convection flow of viscous incompressible fluid through a horizontal channel embedded in non-Darcy porous medium with Brinkman-extended Darcy drag, assuming that temperature of the lower wall varies sinusoidally in the direction of the flow and that of the upper wall is uniform. The governing equations of the flow in terms of vorticity equation and energy equations are simulated using the lattice Boltzmann method together with the finite difference successive over relaxation method. The results are presented in terms of streamlines and isotherms showing the effect of Darcy drag and the Forchheimer drag as well as average Nusselt number. The observation from the present investigation is that the average Nusselt number decreases due to increase of the Darcy and Forchheimer drags, which leads to disappearance of the separated flow that developed for the flow of pure fluid as well as reduces the temperature along the heated region of the lower surface, although the periodicity of the wave propagation remains the same but amplitude of oscillation diminishes.


INTRODUCTION
Considerable interest has been shown in studying of flow of viscous fluid through a porous medium by a significant number of scientists because of its natural occurrence and importance in both geophysical and engineering environments. Another area of nuclear engineering applications in the designing pebble bed reactor which requires a proper understanding of forced convection through packed beds under normal operating conditions and of free convection. Mixed and natural convection flow within the fluid-saturated porous medium is also of interest in a range of other applications; such as water movement in geothermal reservoirs, underground spread of waste, nuclear waste repository and insulation engineering. With this understanding, Koh and Colony (1974) and Wang and Tien (1984), first, analyzed the enhancement of heat transfer by insertion of solid matrices flow through a channel. Enhancement of cooling effectiveness by porous materials in coolant passages investigated experimentally by Koh and Stevens (1975) for low permeability.
Based on the modified Darcy model for transport of momentum Kaviany (1985) studied laminar convective flow through a porous channel surface temperature of which maintained at uniform temperature. Later, Nakayama et al. (1988) considering the Brinkman-extended Darcy model, considering the case of constant heat flux from the walls, studied the effect of the boundary viscous frictional drag and heat transfer characteristics. Combined natural and forced convective flows through a horizontal porous channel connecting two reservoirs had been investigated by Haajizadeh and Tien (1984) and found that small rate of through-flow contributes significant effect on heat transfer across the channel.
The heat transfer characteristics from forced pulsating flow in a channel with uniformly heated surfaces filled with fluid-saturated porous media studied numerically by Kim et al. (1994) considering the Brinkman-Forchheimer-extended Darcy model. Before this analysis, Lauriat and Prasad (1989) investigated the natural convection flow of viscous fluid in a differentially heated vertical cavity filled with porous media giving importance to Darcy-Brinkman- Fig. 1: Flow geometry and coordinates system Forchheimer drag. On the other hand, Ettefagh and Vafai (1988) and Ettefagh et al. (1991) investigated the natural convection flow of viscous incompressible fluid in open-ended cavity with porous obstructing media. Latter, the classical Darcy formulation, the Darcy-Brinkman model, Darcy-Forchheimer model and the Darcy-Brinkman-Forchheimer model has thoroughly been reviewed by Nield and Bejan (2006).
Recently, Hossain et al. (2013) numerically investigated the combined effect of conductionconvection-radiation on the natural convection flow of viscous incompressible and optically thick grey fluid considering the Rosseland diffusion approximation in a porous media with Darcy-Brinkman-Forchheimer drag in a square cavity. Later, Hossain et al. (2013), also investigated the surface radiation effect on the natural convection flow in a fluid-saturated porous medium in a rectangular enclosure using the Darcy-Brinkman-Forchheimer model, considering an exponentiallyvarying surface temperature on the left vertical wall and cool right vertical and top walls. In this study, the bottom wall was considered to be a diffusive gray radiator with emissivity.
In this study we aim to investigate the twodimensional mixed convection of a viscous incompressible fluid flowing through a horizontal channel with spatially periodic temperature boundary conditions at the lower wall and the temperature of the upper wall maintained at ambient temperature, T 0 . The medium is considered to be porous due to Darcy-Brinkman-Forchheimer drag. Wang et al. (1991), Tangborn (1992) and Zhang and Tangborn (1994) carried out their numerical simulations on the unsteady two-dimensional mixed convection flow of viscous incompressible fluid in the absence of Darcy-Brinkman-Forchheimer drag through a horizontal channel with spatially periodic heating from the lower surface. Wang et al. (1991) found a region in (Re, Ra) space where the flow is unsteady and extends as low as Re = 2.0. The above analysis also showed that if the upper wall is insulated, the flow is stabilized and will not become unsteady unless the Reynolds number is raised to Re = 100. These results showed that the flow becomes unsteady when there is a balance between buoyancy and pressure forces. The critical Reynolds number will be lower if the fluid is unstably stratified.
For the present investigation the dimensionless momentum equations that govern the flow of viscous incompressible fluid through a non-Darcy porous media with Darcy-Brinkman-Forchheimer model are transformed to vorticity equations. Lattice Boltzmann Method (LBM) is being applied to simulate the vorticity equation and the related energy equations; whereas, the stream function equation has been simulated applying the Successive Over Relaxation (SOR) finite difference method. The solutions thus obtained are discussed in terms of streamlines, isotherms and average Nusselt number. Further, in this investigation, fluid considered here is water at 20°C for which value of Prandtl number Pr = 7.0. Finally, we have observed the effect of Darcy number, γ and Forchheimer number, Γ, on the formation of separated flows and the temperature distributions in the flow field. We also observed the presence of these numbers on average Nusselt number.

Formulation of physical model:
Consider an unsteady mixed natural convection flow through a horizontal channel of height H filled with fluid-saturated porous medium as shown in Fig. 1. Here we assume that: • All the properties except the density are constant • The density variation is incorporated by invoking the Bousenesq approximation for buoyancy driven flows • The top wall is cold and the bottom wall is maintained at temperature T w (x), that is varies sinusoidal in X, that measures the distance from the entry point down to the right exit point.
Considering the above assumptions, for an isotropic, homogeneous, fluid-saturated porous medium, the governing equations which include the Forchheimer and Brinkman modifications (Lauriat and Prasad, 1989;Hossain et al., 2013;Ettefagh and Vafai, 1988;Hossain and Wilson, 2002;Kaviany, 1985), the energy equation and the equation of continuity may be written as: where, where u, v are the fluid velocity components along the x-and y-axes which are parallel and normal to the channel, respectively. T, ρ, ν, Cp and β are the fluid temperature, density, kinematic viscosity, specific heat at constant pressure and coefficient of thermal expansion, respectively. Further, g is the gravitational acceleration, F, is the empirical constant in the secondorder resistance, κ is the coefficient of thermal conductivity and t is the time.
In Eq. (1) and (2), K is the measure of permeability of the porous medium (a packed bed of spheres), defined by: Here, V f is the volume of the fluid and V C is the control volume. The initial and boundary conditions are given by: Now taking H as the reference length and T 0 is the temperature of the top wall of the fluid in undisturbed channel, as the reference temperature, the dimensionless dependent and independent variables are defined as: , Introducing the variables given in (6) into the set of Eq. (1) to (5) take the following form:  (7) and (8) γ and Γ are, respectively, the inverse of Darcy number and the Forchheimer are defined by γ = 1/(Da Re) and Γ = Fr /Da 1/2 , respectively. It's worthy to mention that these pure numbers, with brevity, as the Darcy drag and Forchheimer drag. Corresponding boundary conditions become: 0 at 0, 2, and 0 at 2 Here, we can see that in the absence of Darcy-Forchheimer drags (that is, for γ = Γ = 0) the present model represents the one investigated by Zhang and Tangborn (1994) and Nakayama et al. (1988). In their analysis Zhang and Tangborn (1994) first expanded the dependent functions by Fourier series in X involving the coupled system of three (including the energy equation) second-order and then the resulting equations solved by using the Chebyshev-Tau method. Time stepping was conducted using Crank-Nicolson for the viscous (linear) terms and Adams-Bashforth for the convective (nonlinear) terms. In this, the twodimensional mesh used for most of the simulations was 32×32 (X and Y directions), though a mesh of 64×64 points has also been used in some cases to confirm the spatial resolution. An example of this simulated result is shown in Fig. 2.
For convenience, we now, eliminate the pressure gradients between Eq. (7) and (8) to find: where, Ω is the vorticity function, which is defined as given below: The energy equation is: Now define: where, ψ is the stream function that satisfies the equation of continuity. Using these in (13) we have: Below we discuss the method applied to finding the solutions of the model presented through Eq. (12)-(16) that satisfying the boundary conditions: In the present simulation we have taken h = 2.

METHODS OF SOLUTION BY LATTICE BOLTZMANN METHOD AND FDM
The Lattice Boltzmann Method (LBM) is a computational and modelling method that differs from traditional numerical methods. It has unique features that other numerical methods do not have due to its micro-particle characterization. Details of this method are available in the book by Mahamad (2011).
Recently, the lattice Boltzmann method has been applied by Alamyane and Mohamad (2010) to investigate the heat transfer enhancement in a channel with extended surfaces mounted above the bottom plate. An investigation has also been made on steady flow and heat transfer in incompressible through a twodimensional constricted channel by Gokaltun and Dulikravich (2010). Shokouhmand et al. (2009) explored the laminar flow and convective heat transfer between two parallel plates of a conduit, fully and partially filled with a porous medium. Viscous and inertial flow resistance effects of a porous medium were also incorporated in the form of force terms in the Boltzmann equation; and a natural convection flow in an open-ended square cavity packed with a porous medium has been simulated by Haghshenas et al. (2010).
The D2Q9 lattice Boltzmann model was used to simulate fluid flow in a two-dimensional channel with uniform grid size of ∆x by ∆y. The lattice Boltzmann equation (known as the LBGK equation) with two relaxation times, Eq. (12) and (14), can be expressed as: And: where, Fr is the source function defined in (28), f k (X,τ) and g k (X, τ) (for k = 0, …., 8) are the distribution function for the vorticity, Ω and the temperature distribution, θ, at time τ in the space grid points X along the direction e i . f k eq (X, τ) and g k eq (X, τ) are the corresponding local dynamic distributions for Ω and θ, c k = ∆X/∆t is the particle speed and ω Ω = 1/τ Ω , ω θ = 1/τ θ ; τ Ω and τ θ are the relaxation times at which the particle distribution functions reache the equilibrium state, which are: Together with: The vorticity and the fluid temperature are determined from the following relations: The evolution of the LBGK model is a relaxation process, through the micro-particle distribution function f i and g i and the equilibrium state f i eq and g i eq to speed up relaxedly, allowing the system to quickly evolve to meet the objective laws of the state. In the LBGK model, the most common is the DdQq series model, where d is the space dimension, q represents the number of discrete speeds.
Here, the local dynamic equilibrium distribution functions of the model D2Q9 are taken into consideration: And: Here C k are unit vectors along the lattice streaming and directions w i the weighted function: And: Figure 3 illustrates the streaming direction and values of w i and C k . The streaming process is as follows: Now from Eq. (22) and (23) In Eq. (24) Φ stands for the vorticity function Ω and the temperature function θ.
The present numerical investigation considers the classical geometry of a rectangular horizontal channel temperature of the upper wall of which is maintained at initial fluid temperature, whilst the lower wall is maintained at a temperature that is sinusoidal with X, that measures the length of the channel from the entry point at X = 0. The fluid enters the channel through the left end, having temperature that is parabolic about the centre line of the channel; the right side is open, the condition at which is assumed to be adiabatic.
Below, the appropriate boundary conditions for the vorticity and the energy function in terms of the distribution function are given: where, θ w = sin(mX), Ω s =-(4u(i,1) -u (i, 2)/ (2dY)). For Eq. (17) we now need to find Fr from: Has been used to find the stream function for certain known vorticity distributions with Ω that assign the optimum value: In Eq. (30) k represents the iteration number. In the present computation, the maximum error is allowed to be 1.0×10 -6 for the SOR. The boundary condition shown in Fig. 3 for ψ = 0 on the solid walls comes from the fact that there is no net flow across these boundaries; at the exit point ߲ψ / ߲Xψ must also vanish along the Y-axis in order to make the fluid motion to the right the mirror image of that to its left. Thus stream function ψ is calculated based on the vorticity distribution by solving Eq. (11) using the Successive Over Relaxation (SOR) method from which we find U and V from the relations given in (15). In the present simulation, the value of Lx has been chosen as 2π. The grid selected for this simulation is 220×80. The present simulated results are compared with those obtained by Zhang and Tangborn (1994) in terms of streamlines and isotherms, taking Re = 5.0, Ra = 20,000 and Pr = 1 in Fig. 4, in the absence of Darc-Frochheimer drag (i.e., γ = Γ = 0). One can see that the present result qualitatively agreed with that of Zhang and Tangborn (1994).
As in Hossain et al. (2013) and Hossain and Wilson (2002), a grid dependence study has been carried out for the present unsteady mixed convection flow through a channel embedded in saturated porous media with Darcy-Forchheimar drag, taking Re = 10, Pr = 1.0, Ra = 20,000, but for γ = Γ = 0 with mesh-points125×40, 185×60 and 220×70. The numerical values obtained for average Nusselt number, Nu av , against time, τ, are shown graphically in Fig. 5. The comparison shows that the choice of meshes 220×70 is better than all other choices, even better than that 51×51 selected by Zhang and Tangborn (1994). For better computational results, the 220×70 mesh has been used throughout the present simulations with maximum time τ = 10 in order to find a steady state solution, as demonstrated in Fig. 6.

RESULTS AND DISCUSSION
In the presented analysis we have considered the effect of Darcy-Forcheimar drag on the unsteady mixed convection flow of a Newtonian fluid through a horizontal channel placed in a fluid saturated porous medium. The temperature of the upper wall was assumed to be constant and that of the bottom surface varied sinusoidally in the horizontal direction. The governing momentum equations were transferred to the vorticity equations. Finally, the vorticity equations and the energy equations were simulated using the Latice Boltzmann Method (LBM) and the stream function equation by the SOR method. The results thus obtained are presented in terms of streamlines, isotherms and heat transfer rate of the walls for different values of the governing parameters. Figure 4a shows the steady streamlines (a) and isotherms (b) for the case Ra = 20000, Re = 5 and Pr = 1 for pure fluid obtained by Zhang and Tangborn (1994). Here the flow looks like a distorted Benard convection and the vertical velocity is quite large. The isotherms show that a plume, which is pushed very slightly downstream, has formed above the hot section of the lower wall. This flow is clearly dominated by the buoyancy force. Almost the same result is shown for the present solution in Fig. 4b. Undoubtedly, from comparison of these two figures, we can claim that the present results agree qualitatively with those obtained by the above authors.
The rest of the results were obtained for Re = 10.0, Ra = 8.0×10 4 and Pr = 7.0 and show the effect of the Darcy and Forcheimer drags in the flow field.
Effect of Darcy-Forchheimer drags on streamlines and isotherms: Steady streamlines and isotherms are shown in Fig. 7a to 7c for values of γ equal to 0.0, 5.0 and 10.0 while the value of Forchheimer drag is absent (Γ = 0). Figure 4a represents the flow for pure fluid for which Pr = 7.0. The results represent streamlines and isotherms, the natures of which are very close to those obtained and shown in Fig. 5b. The difference is only in terms of the intensity of the developed primary (-ve) and secondary (+ve) vorticities. From these figures the numerical values of the rate of intensity for primary vortices are -1.0, -0.42 and -0.23, respectively for values of Darcy drag γ equal to 0.0, 5.0 and 10.0. Further, we can see from these figures that the size of the vorticity reduces with the increase in Darcy drag, as expected, since increase in the Darcy force opposes the velocity of the fluid flow. Now we look at the effect of increase in the Darcy force on the isotherms. As in Zhang and Tangborn (1994), the isotherms show that a plume is pushed very slightly in the downstream regime and forms above the hot section of the lower wall. It can also be seen that the size of the developed plume decreases with the increase of Darcy drag; and it is possible since increased porosity reduces the rate of heat transfer from the hotter regime near the lower plate and also due to decrease in the buoyancy force.  Fig. 6a to 6c, respectively. As in the case of Darcy drag, here we see that increase in the Forchheimer drag reduces the flow rate and the size of the primary as well as the secondary vortices decreases and this results in a decrease in the plume size. The reason is the same as that given for the case of Darcy drag. In both cases, further increase in Darcy and Forchheimer drags will reduce the formation of separated flow in the flow domain. average Nusselt number that in the heat transfer from the lower plate. This is the reason for the nature of the changes in the streamlines and isotherms that what are previously observed.

Effect of Darcy-Forchheimer drags on longer on local shear stress and Nusselt number:
Now we present the numerical values of local shear stressτX and local Nusselt number, Nux, obtained from the following relations: are depicted, respectively, in Fig. 9a and 9b  Effect of Darcy-Forchheimer drags on longer wavelength instabilities: Zhang and Tangborn (1994) also found in carrying out simulations in a periodic domain there is always the danger of creating a nonphysical flow by over-truncating the domain in the periodic direction. The most unstable wavelength could then be longer than the periodicity of the geometry and this instability would not be observed. The simulations in that paper have been tested by these authors and it was found that this type of truncation by doubling the domain size while keeping the boundary conditions fixed. They also investigated the occurrence of this behaviour by decreasing the length of the periodicity of the boundary condition, i.e., by increasing the parameter m. Figure 11a and 11b represent the steady streamlines and isotherms, obtained by the aforementioned authors for the case m = 4, Ra = 15000, Re = 10 and Pr = 1.0.
In this case they found that the flow periodicity matches the boundary conditions that can be seen from the four vortices. As expected, the oscillations occur at higher values of Rayleigh number, Ra and Reynolds number, Re and the periodicity of the oscillating flow also follows by change. Figure 12 depicts the streamlines and isotherms at steady state for m = 4 in absence of Darcy and Forchheimer drag in the flow field which are similar to that of Zhang and Tangborn (1994), i.e., the flow periodicity matches the boundary condition periodicity from the appearance of the four vortices. Difference is only in the size of the and it is because of larger value of Re and Ra. Similar periodicity also occurred in case of the isotherms. From Fig. 9 and 12, depicts, respectively, the streamlines and isotherms for the case when the values of Darcy drag parameter γ equals 10.0 and that of Forchheimer drag parameter Γ equals 10 from which one can see that the periodicity of the oscillation of both the streamlines is the same as in the case of the flow of a pure fluid (Fig. 13).

CONCLUSION
Direct numerical simulations of the unsteady steady two-dimensional mixed convection flow in a horizontal channel immerged in fluid over saturated porous material with spatially periodic heating were carried out using 220×70 mesh points. Various methods were applied in this study, including the Lattice Boltzmann method for the vorticity equation, the energy equation and the successive over relaxation method for the stream function equation. In this study, the base fluid used was water for which the Prandtle number was Pr = 7.0, Reynolds number, Re = 10.0 and Rayleigh number Ra = 8.0×10 6 . Special attention has also be given to the effects of Darcy and Forchheimer drags.
It may further be concluded that, due to the presence of Darcy and Forchheimer drags, the average Nusselt number decreases and this leads to disappearance of the separated flow that happened in case of the flow of pure fluid; this also cause to reduce the temperature of the fluid in the heated region of the lower surface. Moreover, periodicity of the wave propagation remains the same even in Darcy-Forchheimer porous media that happened to increase in the case of pure fluid. Amplitude of oscillation as well as the waviness diminishes in case of higher porosity of the media.