Heat Transfer in Immicible Fluids through a Channel with Porous Beds Bounded by Di ff erentially Heated Plates Using Galerkin ’ s Finite Element Method

Received: November 29, 2010 Accepted: December 10, 2010 doi:10.5539/jmr.v3n2p74 Abstract We analyze the heat transfer in the flow of two viscous incompressible immiscible fluids in a channel with porous beds bounded by differentially heated rigid plates by using Galerkin’s finite element method. Solutions of the governing equations have been obtained by dividing the flow region into four zones applying appropriate matching conditions. The velocity, temperature and the shear stresses, Nusselt numbers are evaluated using finite element analysis and their behavior is discussed for variations in the governing parameters.


Introduction
A system consisting partially of a fluid saturated porous material is known as a composite system, and the convection flow and heat transfer in composite systems form an important class of problems.The convection phenomenon in these systems is usually affected by the interaction of the temperature and flow fields in the porous spaces and the open spaces.The importance of this class of problems is justified both in fundamental and in practical sense with reference to practical thermal engineering applications which stand to benefit, if a better understanding of heat and the fluid flow processes in composite systems is acquired.The following examples are cited; fibrous and granular insulation which occupies only a part of the space between a hot and a cold boundary, fault zones in geothermal systems, cooling of stored grain etc.The heat transfer in flow past a permeable bed was investigated by (Vidyanidhi, V. et al, 1970), (Vidyanidhi, V. et al, 1977), (Rudraiah, N. et al,1977) and (Rudraiah, N. et al, 1976), taking the dissipation effect in the energy equation.Their configuration consisted of a fluid bounded below by a permeable bed.The rate of heat transfer between the fluid and the permeable surface was discussed, making use of slip temperature condition at the interface.(Bhargava, S.K. et al, 1989) investigated the heat transfer in the generalized couette flow of two viscous, inviscid, immiscible fluids bounded by permeable bed making use of Brinkman model (Brinkman, H.C., 1947).Analytical solutions were obtained for the velocity and temperature distribution by dividing the flow region in four zones and applying the appropriate matching conditions.

Formulation
Consider the heat transfer in flow of two immiscible fluids in a horizontal channel of height 2h, bounded by permeable beds of different permeabilities.The outer impermeable walls are maintained at constant temperature.The velocity and temperature field in the porous and clean regions of the channel are assumed to be fully developed.There is a coupling between temperature and velocity in all the zones.The flow takes place under the influence of a constant pressure gradient on each fluid.The entire flow configuration is divided into four zones.The flow geometry described in Figure 1.Zone -I (−h ≤ y ≤ 0) corresponds to a region of clean fluid with viscosity μ 1 bounded below by the porous bed, Zone -III (− (s + h) ≤ y ≤ h) with permeability k 1 , y = −h is the nominal surface separating the clean fluid and porous bed and s is the thickness of permeable material.Zone -II (0 ≤ y ≤ h) contains the flow of clean fluid with viscosity μ 2 (< μ 1 ) , bounded above by the porous bed Zone -IV (h ≤ y ≤ s + h) whose permeability is k 2 , and is separated by the nominal surface y = h.
In the absence of an extraneous force, the following are the equations governing each zone and respective conditions including the interface conditions.

Zone -I :
dV Where V 1 , V 2 , V p 1 , V p 2 are velocities in Zones I, II, III and IV respectively.T 1 , T 2 , T p 1 , T p 2 are the corresponding temperatures.K 1 and K 2 are thermal conductivities of the lower and upper clean fluids respectively.
The interfacial continuity conditions are The boundary condition on the impermeable wall is given by We define the following nondimensional variables where U is the interface velocity, and ρ 1 , ρ 2 are the densities of the lower and upper immiscible fluids respectively.
Making use of these nondimensional variables, the nondimensional equations corresponding to each zone are (dropping asterisks) Zone -I : where are Darcy parameters numbers in zones III & IV respectively is the Brinkman number, and The interfacial conditions and boundary conditions in nondimensional form are The boundary conditions on the impermeable wall in nondimensional form are

Shear Stress and Nusselt Number
The Shear stress on impermeable boundaries in nondimensional form is The rate of heat transfer (Nusselt number) on impermeable boundaries in nondimensional form are

Finite Element Analysis
According to (Zienkiewicz, O. et al, 2005) and (Reddy, J.N., 2005), the underlying mathematical basis of the finite element method first lies with the classical Rayleigh-Ritz and variational calculus procedures.These theories explained why the finite element method worked well for the class of problems in which variational statements could be obtained (e.g., linear diffusion type problems).However, as interest grew in the application of the finite element method to more types of problems, the use of classical theory to describe such problems became limited and could not be applied, e.g., for fluid-related problems.Extension of the mathematical basis to non-linear problems was achieved through the method of weighted residuals (MWR), originally conceived by Galerkin in the early 20th century.The MWR was found to provide the ideal theoretical basis for a much wider variety of problems as opposed to the Rayleigh-Ritz method.Basically, the method requires the governing differential equation to be multiplied by a set of predetermined weights and the resulting product integrated over space; this integral is required to vanish.Technically, Galerkin's method is a part of the general MWR procedure, since various types of weights can be utilized; in the case of Galerkin's method, the weights are chosen to be the same as the functions used to define the unknown variables.The finite element method now employs Galerkin's method to establish the approximations to the governing equations We now apply the Galerkin finite element method to the differential equations ( 12) to ( 19), while dividing the flow region into line elements.The following procedure leads to the finite element equations for velocity and temperature in each zone.

Zone -I
If V K 1 and θ K 1 are the approximation functions of V 1 and θ 1 in the element e k , we define the errors (residual) Where V K 1 and θ K 1 are linear combinations of Lagrange polynomials in terms of the respective local nodal values Where Ψ K j s are shape functions which are given in appendix I.These errors (23-24) are orthogonal to the weight function over the domain e K .Under the Galerkin method, we choose the approximation functions as the weight functions.Multiplying both sides of ( 12) and ( 13) by the weight functions and integrating over the domain e K we obtain Integrating by parts, these line integrals are 25) and ( 26) respectively, we get where We make use of quadratic polynomial approximations, and divide each zone into m i (i = 1, 2, 3) quadratic elements.For computational purposes, we choose two quadratic elements in each zone to assemble the corresponding stiffness matrices in each zone, making use of inter-element continuity, boundary conditions and equilibrium conditions of the secondary variables to obtain the 17 × 17 global matrix in terms of the global nodal values U i (1 ≤ i ≤ 16).Details of the stiffness matrices for each zone are omitted herein for brevity.The corresponding global matrix with respect to the velocity and temperature are discussed below.
Making use of the boundary conditions for the velocity in terms of the nodal values, we obtain The interfacial velocity conditions in terms of the nodal values are given by Equilibrium conditions of the secondary variables are utilized to assemble the matrix equations for the four zones.We obtain the global nodes U i (i = 2...16) , and reduces to a 17 × 17 matrix equation.This 17 × 17 global matrix equation with respect to the velocity can be partitioned in the form where Reynold's numbers in Zone -I and II respectively.
where D 1 , D 2 are the Darcy parameters in Zones III & IV respectively.Solving the matrix equation ( 29), we obtain the solution for U i (i = 2, ...16).
The finite element solutions for the velocity with respect to the four zones are Ψ 1 1 etc. are shape functions under a quadratic polynomial approximation, the details of which are provided in appendix 1.The boundary conditions for the temperature in terms of the nodal values are The interfacial temperature conditions in terms of the nodal values are Equilibrium conditions of the secondary variables at the nodes Making use of the aforementioned conditions and proceeding in a similar fashion, the 17 × 17 matrix equations in terms of nodal values of temperature θ i (i = 2, ..., 16) and where Δ 1 θ , Δ 2 θ , F 1 θ , F 2 θ are column matrices given by G 21 , G 12 and G 22 are similar 8 × 9 and 9 × 9 matrices respectively where Br 1 , Br 2 are Brinkman numbers in Zone-III & IV respectively.Solving the matrix equation ( 31) we obtain the solution for θ i (i = 2, ..., 16) The finite element solutions for the temperature of the four zones Ψ 1 1 etc. are shape functions under the quadratic polynomial approximation, the details of which are given in appendix I.

Discussion
The velocity, temperature are computationally evaluated and the respective profiles in each zone are plotted in Figures ( 2)-( 53) for variations in the governing parameters D 1 , D 2 , R 1 , R 2 , Br 1 and Br 2 .The stresses on the impermeable boundaries and the corresponding rate of heat transfer (Nusselt Number) have been evaluated and summarized in Tables ( 1)-( 4).It is interesting to note that the behavior of the velocity and temperature in case of their porous beds has features different from the case of thick porous beds.
The behavior of velocity, with reference to the variation in permeability of the upper and lower porous beds, may be observed from Figures ( 10) -( 29).When the permeability of upper porous bed corresponds to a Darcy parameter D 2 < 10 4 , any increase in the Darcy parameter D 1 in the lower porous bed reduces the velocity in entire fluid region.However, for D 2 > 10 4 , any increase in D 1 reduces the velocity in lower half of the fluid region while enhancing the upper half of the fluid region.Likewise, by maintaining the permeability of the lower porous bed at D 1 ∼ 10 4 , we see that a lower permeability in upper porous bed creates a higher velocity in the entire fluid region.This is so except at the level y = −1.05,where the velocity transitions from positive to negative.We also observe that by increasing the Darcy parameter in the upper porous bed to D 2 > 2 × 10 4 , the velocity in the fluid region slightly depreciates at all corresponding points while remaining invariant in the upper porous bed.The velocity in all four zones is augmented with an increase in the Reynolds numbers relative to both of the immiscible regions, as seen in Figures (10)-(25).We also observe that, in general, the velocity in the permeable beds reduces with a decrease in the permeability of the beds.This is depicted in [22][23][24][25].
The temperature properties are shown in Figures ( 26)-( 53) for variations in the governing parameters.We find that the temperature increases with increasing D 1 but decreases with increasing D 2 in the clean fluid zones (Figures 26 -44).In the porous beds, with reference to D 1 , the behavior in lower porous bed contrasts with that in upper porous bed.Also the thicknesses of the beds play a significant role on the temperature (Figures 43-53).We may also note that the influence of the temperature with reference to D 2 on the upper porous bed is similar to that of D 2 on the lower porous bed.The temperature in clear fluids as well as porous beds increases with increasing Br 1 and Br 2 , as seen in Figures (32,46,47,40).
The stresses corresponding to variations in the Darcy parameters and Reynolds number are evaluated for different thicknesses of the porous beds (See Tables 1 and 2).We observe that the stresses on the lower and upper plates diminish if the permeability of the porous beds bounded by these plates decrease.The stresses increase with increasing Reynolds numbers corresponding to either of the immiscible fluids.
The rate of heat transfer on the boundary plates for variations in D 1 , D 2 , R 1 , R 2 , Br 1 and Br 2 are presented in Tables 3 and  4. The Brinkman numbers Br 1 and Br 2 indicate the extent to which the viscous heating is important relative to the heat flow resulting from the imposed temperature difference.We find that the Nusselt numbers Nu 1 and Nu 2 depend on the thickness of the porous bed.It is interesting to note that in case of thick porous beds, Nu 1 and Nu 2 increase with increasing D 1 but decrease with increasing D 2 .However in case of thin porous beds, Nu 1 and Nu 2 decrease with increasing D 1 or D 2 .Nu 1 or Nu 2 decrease with increasing Reynolds number R 1 or R 2 on both boundaries, and increases with increasing Brinkman numbers Br 1 or Br 2 ..40 −12.17 −13.78 −8.97 −10.57 −10.41 −11.41 −10.89 −12.50 −15.54 −14.89 0.50 −11.07 −10.68 −4.12 −7.60 −3.41 −2.71 −15.96 −19.44 −22.20 28) are the local (stiffness) matrices of the Zone -I.We can implement the same procedure on Zones -II, III & IV.

Figure 1 .
Figure 1.Schematic diagram of the parallel plate channel flow bounded by porous beds with heat transfer