MICROSCOPIC MODELING OF TURBULENT FLOW OVER AND WITHIN A POROUS BED

P. Prinos¹ and D. Sofialidis²

¹Hydraulics Lab., Civil Eng. Dept., Aristotle Univ. of Thessaloniki, 54006 Thessaloniki, Greece

Tel. & Fax: +30 31 995689, E-mail: prinosp@civil.auth.gr

²Enervac-Flutec Ltd., 57001 Thermi Thessaloniki, Greece

Tel.: +30 31 465013, Fax: +30 31 465056, E-mail: dimitris@flutec.gr 

Abstract: In this study the characteristics of turbulent flow in an open channel with a porous bed are studied numerically. The Reynolds-Averaged Navier-Stokes equations are solved numerically in conjunction with the low-Reynolds k- turbulence model of Launder–Sharma above and within the porous bed. The latter is represented by a bundle of cylindrical rods (with the rods' axes being normal to the main flow direction) of certain diameter and spacing, resulting in high permeabilities, ranging from 35 to 3300 [mm²] and porosities from 0.7 to 0.9. With such an arrangement the microscopic modelling and simulation of flow within and above the porous layer is possible. Main emphasis is given to the effects of relative porous thickness and Darcy number on the flow characteristics over and within the porous bed. Computed velocities indicate the significant influence of the above-mentioned factors on the flow field. Also, computed discharge indicates the capacity of such channels compared to the respective ones with impermeable bed. Computed velocities and turbulence kinetic energy are compared satisfactorily with available experimental measurements taken at the Hydraulics Laboratory of Aristotle University of Thessaloniki.

 Keywords: turbulent flow, porous medium, darcy number, rods bundle, low-Re turbulent model

1  INTRODUCTION

Flow phenomena and the associated momentum transfer near a porous/fluid interface are encountered in various fields (mechanical engineering, environmental hydraulics, geophysical fluid dynamics, among others). In all problems the knowledge of the flow interaction above and inside the porous medium and the momentum transfer at the interfacial region is very important and is required for the water quality aspects. For example, the hydrodynamics of the Ekmann layer in the benthic boundary layer affect the oxygen flux at the water/sediment interface (Svensson & Rahm, [1]). Also, sediment oxygen demand increases linearly with water velocity above the sediments when the velocities are low (Mackentum & Stefan, [2]).

The effect of a porous medium on the flow above it and near the porous/fluid interface has been studied by several investigators (Beavers & Joseph, [3], Poulikakos & Kazmierczak, [4], Rudraiah [5], Vafai & Thiyagaraja, [6], Sahraoui & Kaviany, [7], Alberto Ochoa–Tapia & Whitaker, [8] & [9], Choi & Waller [10]) for relatively low Reynolds numbers (laminar flow). The main problem is the determination of proper interfacial conditions since the porous/fluid interface cannot be considered as a ²no-slip² interface (Beavers & Joseph, [3]).

Studies of turbulent flow near a fluid/porous interface are rather limited, due to the additional difficulties caused by turbulence effects. Studies in MIT (Munoz-Goma & Gelhar, [11], Ruff & Gelhar, [12] and Chu & Gelhar, [13]) have shown that the logarithmic law of the wall is valid for turbulent flow over a porous bed, but the von–Karman constant, , has to be reduced from 0.41 (flow over a smooth impermeable bed) to 0.26. Nezu [14] conducted some experiments and found that  decreases with increasing permeability. Zippe & Graf [15], based also on experimental findings, concluded that the boundary resistance of the tested permeable surface is higher than that of a non-permeable boundary with identical roughness.

Mendoza & Zhou [16] and Zhou & Mendoza [17] presented analytical results for the turbulent flow characteristics over a porous bed and for the velocity distribution within the porous bed, respectively. Also, Li & Garga [18] presented analytical results for the turbulent seepage flow occurring at the transition zone (between the fluid zone and the pressure seepage zone) of gravel river reaches or non-conventional rockfill spillways. However, the whole analysis was based on known velocity and shear stress at the interface (top of the transition seepage zone) from the main channel flow characteristics (Li, [19]).

In the present study the characteristics of turbulent flow in a two–dimensional open channel with porous bed, are studied numerically. The Reynolds-Averaged Navier-Stokes equations (RANS) are solved in conjunction with a low-Reynolds k- turbulence model (Launder & Sharma, [20]) above and within the rods' bundle, as the flow may exhibit laminar regions inside the porous medium. The equations are solved with a finite-volume method using the FLUENT 5.3 code [21]. The porous layer is actually a bundle of cylindrical rods (their axes being normal to the flow direction) of certain diameter and spacing resulting in porosity which ranges from 0.7 to 0.9. Emphasis is given in the effects of relative porous thickness and the Darcy number, Da(=K/h΄², where K=permeability and h΄=flow depth over the porous layer) on the flow characteristics over the porous bed. Computed mean velocities and turbulence parameters indicate the significant influence of the above–mentioned factors on the flow. Also, computed discharge over the porous region is compared with the corresponding one of channels with impermeable bed at the location of the fluid/porous interface.

2    GOVERNING EQUATIONS – NUMERICAL PROCEDURE

Using the microscopic approach, the RANS equations are solved in the flow region above the porous medium, as well as inside it. The medium, of given permeability, K, and porosity, , is simulated as a bundle of cylindrical rods (Fig. 1). The geometry in Fig. 1 is repeated infinitely in the flow direction, x. For incompressible, two-dimensional flow the RANS equations are:

Continuity:                                                     (1)

x-momentum:                  (2)

y-momentum:                  (3)

where U, V=components of velocity in the x- and y-direction, respectively; P=effective pressure (the difference between the static and the hydrostatic pressure); , , =Reynolds stresses; ρ=fluid density; and ν=kinematic viscosity. The eddy-viscosity concept is used for modeling the Reynolds stresses appearing in Eqs. (2) and (3). In general:

                        (4)

where ν_t=eddy viscosity; δ_ij=Kronecker delta; and k=turbulence kinetic energy.

The Launder & Sharma [20] low-Reynolds k-ε model is used for calculating ν_t. The use of a low-Reynolds model is necessary in regions where turbulence is damped (porous and near–wall regions). The flow is resolved up to the walls without special treatment (wall functions), with the condition that the first grid point must lie well-inside the viscous sublayer (y⁺<2). Hence, the flow is calculated in both the viscous and turbulent region. The eddy viscosity is related to k and its rate of dissipation, ε, through the Kolmogorov-Prandtl relationship:

ν_t=c_μf_μ(k²/ε)                              (5)

where c_μ=0.09; and f_μ=exp[-3.4/(R_t/50)²]=damping function accounting for low-Re number and wall-proximity effects. The following transport equations for k and ε are solved:

k-equation:           (6)

ε-equation: (7)

where P_k= (¶U_i/¶x_j)=production rate of k due to shear; R_t=k²/(νε)=turbulent Reynolds number; f_ε1=1; f_ε2=1-0.3[exp(-R_t²)]; c_ε1=1.44; c_ε2=1.92; σ_k=1.0; and σ_ε=1.3.

The interaction of turbulent flow above and inside the porous bed may result in turbulence penetration into the upper part of the porous medium, which in turn results in increased turbulence levels (shear stresses and intensities) at the interfacial region (flow/rods' interface). Increased shear is expected to bring a reduction of the velocity in the flow above the rods and a respective increase of velocities in the upper porous area. In other words, turbulence promotes the momentum exchange across the interface, compared to the situation of an impermeable wall located at the interface. Such reduction (from the impermeable wall case) of the velocities above the porous layer is expected to reduce the discharge capacity. Opposite trends have been observed for laminar flow over and within a porous bed. Beavers & Joseph [3] have used the momentum equation for fully developed laminar flow and an empirical fluid/porous interfacial condition and have calculated an increase in mass flow over the porous medium with regard to that with impermeable bed.

Fig. 1     Flow geometry for various porosities

Fig. 2    Numerical domain, and boundary conditions

Fig. 3    Representative mesh used in the computations

The two-dimensional RANS equations in conjunction with the Launder & Sharma [20] low- Reynolds k-ε model are solved with the code FLUENT5.3 [21]. The grid was dense enough near the walls, as indicated by the low y⁺ values (maximum computed value, y⁺=1.05 at the top of the upper rod). Provision was made in order to have at least 2¸5 nodes located inside the viscous sublayer (y⁺<2). The computational domain, as well as the notation, co–ordinate system and boundary conditions are shown in Fig. 2. The geometry is assumed to repeat infinitely in the streamwise direction, x. The numerical mesh, constructed by triangular and quadrilateral elements (around the rods), is shown in Fig. 3. The grid size is given in Table 1.

Periodicity conditions apply in the direction of the mean flow extracted for a given channel slope (given streamwise pressure gradient). At the free surface symmetry conditions are applied since the Froude number is relatively small and there is no surface variation in the flow direction. At the solid boundaries (bottom impermeable wall and rods' surfaces) no-slip boundary conditions are applied.

The main characteristics of the numerical procedure used in FLUENT 5.3 are: (a) Discretization. Second Order for P and QUICK for U, V, k and ε. PISO for velocity–pressure coupling. (b) Number of Iterations. Approximately 50,000 for complete convergence. (c) Convergence. The periodic condition decelerates convergence, hence the large number of iterations. Double precision was used and all normalized residuals were dropped below 10^–8 before convergence was achieved. Convergence was monitored at the periodic and the middle (x=0) planes, in terms of ²frozen² average values of P, U and k.

3    ANALYSIS OF RESULTS

Three test cases with highly permeable beds were studied resulting in porosity, φ, (ratio of volume of fluid to total porous medium volume) of 0.7144, 0.8286 and 0.8929. In all cases three rods were assumed (having diameter, D=10 [mm]) and the height of the porous region, h, was kept constant to 55 [mm]. The geometric and hydrodynamic characteristics are given in Fig. 1 & 2 and Table1. Porous permeability, K, was estimated by a method described in Bird et al. [22]. The channel slope was kept equal to that used in the experiments (S_o=0.002), hence, –dP/dx=ρgS_o=19.62 [N/m³] (where g=9.81 [m/s²]=gravity acceleration).

                  Table1    Geometrical and hydrodynamic data

Fig. 4a, 4b and 4c show the contour plots for normalized U (with U_int, the velocity at the interface, y=0.0), k⁺(=k/U_*²), and turbulent shear stress, ⁺(= /U_*²), respectively, for φ=0.8286 and Da=0.1643 (U_*=(gh΄S_o)^1/2=friction velocity at y=0.0). It is shown that in the porous region close to the impermeable bed the velocities are quite low, especially in the region between the rods, indicating laminar flow. Close to the interface, velocities increase and the low-velocities region in front of and behind the rods has been reduced due to the flow interaction. Mean velocities in the region above the interface are only up to 3.5 times the interfacial velocity indicating the intensive momentum exchange between the fluid and the porous regions. Similar observations are made for the other two quantities, k and . In the lower part of the porous region both quantities are quite low indicating laminar flow, while in the upper part these quantities are increased implying penetration of turbulence.

Figs. 5a & b present the U profile in wall co–ordinates, y⁺(=yU_*/ν) and U⁺(=U/U_*), for relative porous thickness, h/H, equal to 0.4400 and 0.5238, respectively. In the same plot the viscous sub-layer (U⁺=y⁺) and the fully–turbulent (logarithmic) layer (U⁺=(1/κ) lny⁺ + C, where κ=von Karman constant=0.41 and C = integration constant=5.29) for flow over an impermeable wall are also shown. It is evident that the velocities above the porous region are much lower than the respective ones for flow over an impermeable bed. This is due to the action of turbulent shear stress and the penetration of turbulence in the porous region, which reduces the mean velocities above the interface. Near the interface, velocities are shown to deviate from the logarithmic law indicating that the flow near the rods is not of boundary-layer type. Similar trends are observed in the experimental velocity distribution, which although indicates higher velocities than the computed ones, they are still much lower than those over an impermeable bed. The disagreement between experiment and simulation may be due to the limited streamwise length of the porous region in the laboratory channel and therefore to the incomplete development of the flow [23]. It is shown that with a decrease of Da number (decrease of permeability) the velocity above the interface increases and approaches the value of that over an impermeable bed. With a decrease of permeability (decrease of x-spacing of the rods, S) the penetration of turbulence in the porous medium also decreases and therefore the effect of the turbulent shear stresses on the mean flow field becomes weaker. The above effects are more pronounced with increasing relative porous thickness, h/H, i.e. with decreasing flow depth over the porous region the penetrated turbulence to the porous region decreases and hence mean velocities over the porous region are less affected.

Figs. 6a & b show the vertical variation of k⁺ above and within the porous region (the vertical solid lines denote the interface). It is shown that the k distribution is more uniform within the flow depth above the porous region than that observed over an impermeable bed (empirical relationship developed by Nezu and Nakagawa [24]). Computed and experimental k are in satisfactory agreement. Penetration of turbulence to the upper part of the porous region (y/h΄<-0.1) results in increased levels of k, while k is extinguished in the remaining part of the porous layer. Penetration of turbulence into the porous region is shown to be more intensive with increasing Da number (increasing permeability).

Fig. 5    Velocity distribution above the porous region. (a) h/H=0.4400, (b) h/H=0.5238

Fig. 6    Effect of darcy number on k distribution. (a) h/H=0.4400, (b) h/H=0.5238

Fig. 7    Effect of darcy number on distribution. (a) h/H=0.4400, (b) h/H=0.5238

Fig. 7 shows the profile of the turbulent shear stress ⁺ above and within the porous region. It is shown that at the periodic plane (between the two rods) turbulent shear stresses are increased over the whole flow depth above the porous region. In particular, - is higher than U_*², indicating again the mixing layer type of flow rather than the boundary layer type. Again the penetration of turbulence into the upper part of the porous region is indicated by the increased levels of - in the upper part of the porous layer.

Figure 8 shows the effect of Da number (permeability) on the channel discharge indicating the highly reduced discharge capacity of channels with a porous bed. The discharge in channels with highly permeable bed is shown to be in the order of 30% to 45% of (discharge with an impermeable bed).

Fig. 8    Variation of channel discharge with darcy number

4    CONCLUSIONS

The characteristics of flow over and within a highly porous bed have been studied numerically using the Reynolds-averaged Navier-Stokes equations in conjunction with the low-Re k-ε turbulence model of Launder & Sharma [20]. The Da number is varied from 7.3´10^-3 to 3.6752 while the relative porous thickness from 0.4400 to 0.6471. The following conclusions can be derived:

(1) Mean velocities over the porous bed are shown to decrease considerably (with regard to the ones over an impermeable bed) due to the intensive action of turbulence near the fluid/porous interface. Near such an interface the velocities deviate from the law of the wall exhibiting a reduced value of constant C.

(2) Mean velocities over the porous bed decrease with increasing Da number (increasing permeability).

(3) Penetration of turbulence to the upper part of the porous region results in increased levels of turbulence kinetic energy and shear stress in the permeable part below the fluid/porous interface.

(4) The discharge capacity of channels with a highly porous bed is reduced significantly and is shown to be in the order of 30% to 45% of that over an impermeable bed.

References

[1]    Svensson, U., & Rahm, L., (1991) “Towards a Mathematical Model of Oxygen Transfer to and within Bottom Sediments”, J. of Geophysical Res., 96, pp. 2777–2783.

[2]    Mackentum, A. A., & Stefan, H. G., (1998) “Effect of Flow Velocity on Sediment Oxygen Demand: Experiments”, J. of Environmental Engrg., ASCE, 124(3), pp. 222–230.

[3]    Beavers, G. S., & Joseph, D. D., (1967) “Boundary Conditions at a Naturally Permeable Wall”, J. Fluid Mech., 30(1), pp. 197–207.

[4]    Poulikakos, D., & Kazmierczak, M., (1987) “Forced Convection in a Duct Partially Filled with a Porous Material”, J. of Heat Transfer, ASME, 109, pp. 653–662.

[5]    Rudraiah, N., (1985) “Coupled Parallel Flows in a Channel and a Bounding Porous Medium of Finite Thickness”, J. of Fluids Engrg., ASME, 107, pp. 322–329.

[6]    Vafai, K., & Thiyagaraja, R., (1987) “Analysis of Flow and Heat Transfer at the Interface Region of a Porous Medium”, Int. J. of Heat Mass Transfer, 30(7), pp.1391–1405.

[7]    Sahraoui, M., & Kaviany, M., (1992) “Slip and No–Slip Velocity Boundary Conditions at Interface of Porous, Plain Media”, Int. J. Heat Mass Transfer, 35, pp. 927–943.

[8]    Alberto Ochoa-Tapia, J., & Whitaker, S., (1995) “Momentum Transfer at the Boundary between a Porous Medium and a Homogeneous Fluid–I Theoretical Development”, Int. J Heat Mass Transfer, 3(14), pp. 2635–2646.

[9]    Alberto Ochoa-Tapia, J., & Whitaker, S., (1995) “Momentum Transfer at the Boundary between a Porous Medium and a Homogeneous Fluid–II Comparison with Experiment”, Int. J Heat Mass Transfer, 3(14), pp. 2647–2655.

[10]   Choi, C. Y., & Waller, P. M., (1997) “Momentum Transfer Mechanism for Water Flow over Porous Media”, J. of Environmental Engrg., ASCE, 123(8), pp. 792–799.

[11]    Munoz-Goma, R. J., & Gelhar, L. W., (1968) Turbulent Pipe Flow with Rough and Porous Walls, Report No. 109, Hydrodynamics Lab., Dept. of Civil Engrg., Massachusetts Inst. of Technology, Cambridge, Massachusetts, U.S.A.

[12]    Ruff, J. F., & Gelhar, L. W., (1970) Porous Boundary Effects in Turbulent Shear Flow, Report No. 126, Water Resources and Hydrodynamics Lab., Dept. of Civil Engrg., Massachusetts Inst. of Technology, Cambridge, Massachusetts, U.S.A.

[13]    Chu, Y. H., & Gelhar, L. W., (1972) Turbulent Pipe Flow with Granular Permeable Boundaries, Report No. 148, Ralph M. Parsons Lab. for Water Resources and Hydrodynamics, Dept. of Civil Engrg., Massachusetts Inst. of Technology, Cambridge, Massachusetts, U.S.A.

[14]    Nezu, I., (1977) Turbulent Structure in Open–Channel Flows, PhD Thesis, Kyoto Univ., Japan.

[15]    Zippe, H. J., & Graf, W. H., (1983) “Turbulent Boundary–Layer Flow over Permeable and Non–Permeable Rough Surfaces”, J. Hydraulic. Res., IAHR, 21(1), pp. 51–65.

[16]    Mendoza, C., & Zhou, D., (1992) “Effects of Porous Bed on Turbulent Streamflow Above Bed”, J. of Hydraulic Engrg., ASCE, 118, pp. 1222–1239.

[17]    Zhou, D., & Mendoza, C., (1993) “Flow Through Porous Bed of Turbulent Stream”, J. of Engrg. Mechanics, ASCE, 119(2), pp. 365–383.

[18]    Li, B., & Garga, V. K., (1998) “Theoretical Solution for Seepage Flow in Overtopped Rockfill”, J. of Hydraulic Engrg., ASCE, 124(2), pp. 213–217.

[19]     Li, B., (1990) “Characteristics of Flow in Rough Channels with Permeable Bed”, Proc. of the 7^th Congress APD-IAHR, Chinese Association for Hydraulic Research, pp. 1–7.

[20]    Launder, B. E., & Sharma, B. I., (1974) “Application of the energy-dissipation model of turbulence to the calculation of flow near a spinning disk”, Letters in Heat and Mass Transfer, 3, pp.269-289.

[21]   FLUENT (1998) Fluent 5 User’s Guide, Fluent Inc., Lebanon, N.H., USA

[22]    Bird, R. B., Stewart, W. E. & Lightfoot, E. N., (1960) Transport phenomena, Wiley, New York.

[23]    Prinos, P., Sofialidis, D., & Keramaris, E., (2000) “Characteristics of Open Channel Flow Over a Porous Bed”, submitted to J. of Hydraulic Engrg., ASCE.

[24]    Nezu, I., & Nakagawa, H., (1993) Turbulence in open channel flow, IAHR monograph, Balkema Pub.