CFD Analysis for Turbulent Flow within and over a Permeable Bed

Abstact The main objective of this study is to improve our knowledge about the velocity profile and turbulence within and over a permeable bed. The study was using computation fluid dynamics (CFD) methodology to simulate the studied cases. It includes a detail analysis for two-dimensional fully developed turbulent flow over and through a permeable bed. Five different cases were simulated numerically. The analysis is set for the three flow zones (free stream, porous, and interface). The detailed two d imensional flow simulat ions were subsequently validated using previously published results, then it was specially averaged to overcome the heterogeneity of flow. The focus in this study is on the effect of porosity and free stream thickness on longitudinal and vertical velocities in different flow zones. On the basis of this study results, it is shown that the flow velocity within the porous zone increases with bed porosity, and decreases with increasing the water depth. It is also confirmed that the turbulence parameters (turbulent kinetic energy, turbulent dissipation rate, and turbulent kinetic energy production) penetrate practically throughout the whole porous layer to reach maximum values at the interface then decreases smoothly to minimum at the water surface.


Introduction
In general, a permeab le bed has recently treated analogously to an impermeab le bed. The velocity distribution, and flo w resistance coefficients were derived irrespective of bed porosity. However, significant interaction processes occur between the flow above and through the porous zone wh ich depends on the permeability of the surface. These interactions influence a non-zero velocity at the permeable zone, and turbulence exchange of mass and mo mentu m between the two flo w zones. The exchange process is responsible for additional shear stress near the boundary [1].
Nowadays, several studies gave a sight on the velocity distribution within, and close to the porous zone, especially within the interface zone. Such as [2] published that the timeaverage stream-wise velocity increased much more near the free surface than near the permeable bed. It was found experimentally that the velocity profile over the suction zone consists of two parts; an inner boundary layer where the shear stress changes rapidly, and a logarithmic law is applied for the zone above the first layer in which the flo w is unaffected by the presence of suction [3]. The logarith mic zone at the free stream has been widely applied such as [4][5][6][7].
It was found that there is a significant increase in the near bed velocity and a velocity reduction near the water surface, resulting in the format ion of uniform velocity distribution [2].
Although, the vertical velocity distribution is an important variable to be studied and analyzed as it gives a good indication about the fluid movement in the interface region, the previous studies focus on the longitudinal velocity as most of it deals with the free strea m, where no vertical velocity exist.
Fully developed turbulent flow through and over permeab le bed is characterized by the mo mentu m transfer occurring fro m the faster flow at the free stream and the slower flo w in the porous zone. Th is phenomenon occurs at the interface zone, where a shear layer develops. It is developed with the secondary currents generated at the corner regions of the cross section, bed-generated turbulence, and free surface effects. These phenomena result in a very complicated flo w field that studied by lots such as; [5] and [8][9][10]. ϕ δ . 2

*
Cd U U v > < = (1) It was shown in that turbulence can penetrate into the pore space between the obstacles. The interface zone can therefore be defined as the region of the porous mediu m and free stream influenced by the free stream flow. The thickness of the interface zone is called the penetration width.
A dynamic relationship between the roughness layer thickness and (C d α) -1 using data from d ifferent obstructed flows, and found that δ v = 0.333(C d α) -1 [12]. This indicates that the extent to which transport length scale of the thickness (C d α) -1 determines the vert ical transport length scale δ v .
Turbulence is three-dimensional time dependent non-linear phenomenon. The main characteristic of turbulence is the transfer of energy from larger to smaller length scales. Large energy fro m the mean flow is extracted by the large eddies and then transferred to neighbouring eddies of smaller scale. The s mallest eddies then lose energy through the action of viscous dissipation rate, which finally convert into thermal energy.
Turbulence can be described by several variab les, such as, the turbulent kinetic energy, and turbulent dissipation rate. The turbulent kinetic energy is characterized by root-mean -square of velocity fluctuations. Due to the difference between velocities in the free stream and porous zone, the turbulence increases at the interface zone.
Several studies published about turbulent kinetic energy within and over permeab le beds such as [10]. The effect of suction rate on TKE was studied by [15] and found that decreasing the suction rate up to 9% decreases the TKE mo re than four times. [16] found that within the range of Reynolds number of their analysis, the penetration of the turbulence remains independent of the Reynolds number. Further, the existence of an ext ra production term in the TKE-equation is important in the porous region.
Still lots of detailed analysis needed in order to overcome all the factors that affect the turbulent flow within and over permeab le bed such as the effect of flow geometries, and or the effect of the density of the porous zone. Fu rther, literature contains lots of studies about the flo w features and factors affecting them. However, because of the complexity of the system within the free stream and the porous zones, up to the author's knowledge there is a leakage of knowledge at this area of study. Further, few studies in this area used spatially averaged technique to overcome flow heterogeneity due to both turbulence and presence of solid obstacles.
The main objective of this paper is to study the characteristics of the mean flow and turbulence within and over a permeable bed in detail. The whole flow do main can be divided into the three important flo w zones: free stream, porous layer, and the interface between them. The simu lation results include mean velocities and various turbulence quantities (turbulent kinetic energy (TKE), turbulent shear stress, TKE production and dissipation).

Governing Equation
In general, t wo equation models are widely used in present engineering calculations to solve Reynolds Averaged Navier Stokes (RANS) equation. It gives independent scales (eliminates the need for specifying the turbulence length scale) in calcu lating the turbulent viscosity (μ t ). The numerical solution fo r the 2-D model used in this study is K-epsilon. Within the K-epsilon model there are a three solution models available Standard, RNG, and Releasable models. The solution in this study is by two equation model using k-epsilon (k-ε) model wh ich is widely used in literature.
The two equations solved in K-epsilon model are the turbulent kinetic energy, and the turbulent dissipation rate: The dissipation ε: where : In the present study, the commercial CFD software Ansys-Fluent is used. In common with many other codes, Fluent solves the Navier Stokes equation numerically using the fin ite volume method on an unstructured grid. The Navier stokes equations consist of a continuity and three momentum equations, which are based on the conservation of mass and mo mentu m respectively ( [17]).

The Simulated Cases
The study was a numerical simu lation for fu lly developed turbulent flow zone over and through permeable bed. The data for all simulated cases are defined in [18]. They carried out their experimental and nu merical study for three different bed porosities. These cases involved turbulent flow over and within three and four layers of rods. These rods were fixed at different distances to enable to include the porosity effect on turbulence. Reference [18] co mputed the flow structure above and within the porous mediu m using a model based on the Reynolds -averaged Navier-Stokes equations. It was solved using Fluent software with the k-epsilon method. In their simu lation, a periodic boundary condition was used. The periodic boundary conditions may do not accurate results within the porous zone as the flow passes around the rods from vortices the may be cut before it is comp lete and the flow pattern before and after the rod is not the same which is the basic in periodic boundary condition computation. In this study long flu me is built for the studied cases. This enables us to reach fully developed turbulent flow in condit ions as the real experiments.  The geometries of the cases were built in Gambit software. The permeab le bed was consisting of a bundle of rods mounted at certain locations at the bottom of the bed. The water flows over and through these rods. Five different cases were selected with the geometries shown in Table ( 1). Figure  (1) shows the arrangement of rods bundle.
Then it is exported to Ansys12 (Fluent) for simu lation. A sensitivity analysis has been done for spar30 case, as it can be computed with larger mesh size, and as the computation time in Fluent is large, and the expected accuracy will be less for larger mesh. Therefore, the dense cases have not been mesh tested. To test the mesh size affect on the solution. Case spar30 was meshed three times, each with different mesh sizes. The mesh for all cases studied was triangular. Each case having nearly a double mesh size of the previous smaller one, with 0.45 mm for the s mallest size. The testing was for the results that show mesh independence. The results show that the mesh size of 0.98 and 0.45 mm leads to the same results. For this the mesh size of 0.98mm was selected for the spars cases. The exported cases from Gamb it are used in simulation in Ansys12 (Fluent) software, with the flowing boundary condition, and solution controlled methods: The numerical model selected was kε with one of the two models mentioned at section. Under the solution methods the pressure-velocity coupling scheme option selected is "Simp lec" with skewness correction of 2. Simp lec scheme imp roves convergence for pressure-velocity coupling [19]. The spatial discretization selected are; the gradient is "Least Square cell Based", the pressure, mo mentu m, turbulent kinetic energy, and turbulent dissipation rate are all "Second Order Up wind".
Under the solution methods the pressure-velocity coupling scheme option selected is "Simplec" with skewness correction of 2. Simplec scheme improves convergence for pressure -velocity coupling (Fluent user Gu ide, Chapter 18 in www.Fluentusers.com). The spatial discretization selected are; the gradient is "Least Square cell Based", the pressure, mo mentu m, turbulent kinetic energy, and turbulent dissipation rate are all "Second Order Upwind". The appropriate boundary conditions for the p resent problems are set as a uniform flow at the entrance with in let velocity, a no slip condition at the impermeable wall, symmetry at the free stream region, outlet pressure at the exit. The interior edges identified as interior zone in the simulat ion. The simu lation was consumed a long time and large number of iterations to reach fully developed turbulent flow shown in Table 2. The fully turbulent flow is the zone were the flow features do not change with distance. In this study it was tested at each 100 mm distance for mo re than five d ifferent x-locations.

Model Testing and Vali dation
In this study, the values of various flow parameters are kept the similar to that published [18] to assure correctness and consistency of the simulated results. The results simu lated were tested for their validation [18]. Further, t wo simu lated results fro m spars cases were tested with [16] who simu lated these results using macroscopic technique to solve the (RANS) equations. Figure 2 shows a comparison between published and present data for normalized turbulent kinetic energy (TKE + ) vs. normalized bed height for a mid-line between t wo colu mn rods. Fro m the figure, it can be concluded that the simulat ion is with reasonable results that is superior closer to experimental results than previous studies at the interface between the free stream and the porous zone.

Results and Discussion
The work described here was undertaken to provide more detailed information about the velocity field and turbulence within the fully developed turbulent region. The data that exported fro m Fluent was irregular. Further, it was heterogeneous. To overcome these problems, a Matlab program was built for each case to convert the data to regular and to spatially averaged to overcome the flow heterogeneity.
Spatial averaging was for conventional bed-parallel volu mes (centre to centre colu mn rod) for obtain ing vertical properties of velocity co mponents. The spatially averaged quantity of a flow variab le ψ is defined as: Where: h is the flow zone thickness.

The Velocity withi n and over the Permeable Bed
The flo w velocity through and over the porous zone is shown as contours at figure 3. Th is figure is fo r longitudinal velocity for two different porosity cases. From the figure, it could be noticed that the overall flow pattern is similar for both shown cases.
The mean velocities show the expected boundary-layer profile in the free stream zone, and a much slo wer flow in the porous zone. A pipe-like velocity profile fo rms between each two successive rows of cylinders, with a faster and mo re developed flow in case of sparse cylinder packing, co mpared to the dense case. CFD provides details of the flow spatial structure between the cylinders, where the difference between the sparse and the dense packing is most pronounced. As expected bed porosity has a majo r influence on the flow velocities within the porous layer. The longitudinal flow velocit ies across the overall flow profile including both free stream and porous zone are shown in Figure 3. For co mparison, each Figure shows the results for the 'sparse' case (top) and 'dense' case (bottom). The overall flo w pattern for both cases are similar: as expected, velocity magnitudes are much larger in the free stream than within the porous zone; furthermore, within the free stream, the velocity profile has the shape typical for the boundary layer flo ws, whereas with in the porous zone velocity profile is similar to the p ipe flow. The 'p ipes' within the porous zone for 'sparse' and 'dense' case are very different: in the former case they have much bigger diameter, wh ich penetrates into the horizontal gap between the two adjacent rods; in the latter case the 'pipe' is very narrow and it maintains practically constant diameter. The 'sparse' case shows another interesting feature -the maximu m flow velocity in the highest pore is somewhat smaller than in the second pore. This is in agreement with the experimental results of Pokrajac et al. [20], who prov ide a detailed d iscussion of the physical mechanis ms behind this counterintuitive result.  Table 3. Fro m the table it can be concluded that a large discharge difference exists between the porous and free stream zones, even when the area of flow is smaller than that at the porous zone. Fu rther, as expected a co mparison between the studied cases shows that the discharge difference between the free stream, and the porous zone, increases with the density of the rod packing. Furthermore, the larger is the free stream thickness, the larger is the discharge difference. This result is in agreement with experimental results of [21].  In mo re detail, Figure 4, and 5 show the flow vectors coloured by longitudinal velocity for both spar30, and dens30 cases. The flow area over, and underneath the rods are smaller than between two colu mn regions. This decrease in the cross sectional area, increases the velocity. Beyond this region, the cross sectional area increases, so the flow velocity decreases. This pattern obeys fluid dynamics laws "velocity is inversely proportional to flow cross sectional area". At spar30 case, there is sufficient distance for the vortices to form two co mplete vort ices behind the rod and then the flow return back before it reaches the second rod. However, the scenario within the dense cases is different, as due to the small distances between the rods the flow is nearly stagnant in between the rods. As there is no sufficient spaces, the vortices forms at the upper and lower regions of the rods at the centre side between each two-co lu mn rods. In this study, it is noticed that within the interface zone, the velocity is with wave patterns. These waves can be seen at the contour plots in Figure 3. Th is is due to the high flow d ifference between free stream and the upper part of the obstacle, where there the velocity is nearly zero.  between the adjacent rows of the rod. Velocity magnitudes in the pore are smaller than in the second pore for both sparse and dense case. In the former case, the geometry of the two pores is also different, especially in the upper part. As stated earlier, this is the effect of the surface flow wh ich penetrates between rods and changes the geometry of the flow in the first pore.
In closer insight, for the flow features (velocities, and turbulence) were averaged in a Matlab program. The longitudinal velocity within the bed varies from zero close to the solid surface to the highest at the med point between two rows of rods. Figure 7 shows a co mparison between the spare and dense cases for normalized longitudinal velocity distribution within the porous zone. Both cases were normalized (dimensionless) to overcome the geo metry difference between both cases. The velocity was normalized by dividing it by shear velocity (U/U * ). The vertical distance within the porous zone was made dimensionless, using[h* rod = (z + d/2)/h rod -1], where h rod is the rod centre to centre vertical distance in order to align the vert ical position of the rods in the dense and sparse cases. This gives the position of rods 1, 2, 3 as h * rod = -1, -2, -3, respectively. Fro m figure 7, and 8, it can be noticed that as the porosity increases, the bed longitudinal velocity, also increases. Further, for the sparse case, the velocity is small but still positive in the region between the two columns. On the other hand, the dense case shows a stagnant flow close to the rods, that reaches negative values at the upper part of the top rod. The free stream thickness affects the flow within the porous zone porous zone. Figure 9 shows the effect of free stream thickness on velocity within porous zone, a, and b longitudinal velocity for spars, and dens cases respectively, and c is for vertical velocity. Fro m the figure, it can be concluded that the higher the free stream water thickness the lower flow velocity within the porous zone.   be seen that within the porous zone, the velocity changes fro m negative to positive at different parts of the porous bed with a highest value at the interface zone. This feature is due to the change in flow direct ion because of presence of rods. At the interface zone, the flow passing over the upper rod meets with the high speed at the free stream. Both the differences in speed and flow direction between the upper and lower zones of the bed increase the turbulence generation at the interface zone. However, at the free stream the vertical flow decreases to reach zero as the flow changes to a longitudinally uniform flow.  Figure 11 shows the vertical velocity contours, and distribution with the flo w stream for the top rod. Fro m the figure it can be estimated that the spar30 case showed within the porous zone the vertical velocity is fluctuating. It is changing around the rod. It has a high value at the upper inflow, and at the lo wer outflow side of the rod. The high values are due to the change of flow direction because of the presence of the obstacle.
The penetration width (roughness layer th ickness) has been previously studied in several publications such as [22] who suggested that the penetration width is defined as the distance from the surface of the porous zone to the point where shear stress decays to 10% o f its maximu m value. Applying this proposition at the simulated cases resulted penetration thickness values shown in Table 4. Fro m the table it can be concluded that the larger the porosity the lower the penetration width. This gives indication that the transient layer between the free stream (fully developed turbulence), and the porous zone is affected with the free stream thickness. This result is very important for environmental engineering studies at open channels for pollutant migration fro m flo w zone to another.
The dimensionless (δ v /d) is plotted against scale parameter (C d αd) -1 in Figure 12. Where C d is the drag coefficient, α is the frontal area per unit volu me, co mputed for the top rod. Fro m the results it can be concluded that the larger the porosity the larger the penetration width. These results agree with the experimental results of [22][23][24]. Further, the higher the free stream thickness showed a larger penetration layer thickness, which disagree [8].   Figure 13 shows the contour plots and vertical profiles of TKE for sparse and a dense case. For both cases, TKE has a clearly defined maximu m close to the interface. This is due to the difference between velocities in the free stream and porous zone, which creates the most intense shearing at the interface. W ithin the free stream, TKE decreases with distance from the bed. Overall, the TKE pattern agrees with different published data such as; [16] and [18] Turbulence within the porous zone is heterogeneous and characterized with very s mall values of TKE. For the sparse case, TKE within the porous layer steadily increases towards the interface. However, the dense case has higher TKE values close to the interface where TKE reaches maximu m. In this region, the TKE value for the dense case is nearly double that of the sparse case. This is probably due to the higher difference in velocity between the free stream, and the porous zone, wh ich enhances shearing and hence also generation of turbulence. within the porous zone (bottom) for both flow depths the sparse case has higher TKE in both free stream and porous zone compared with the dense case. Further, higher free stream thickness shows a higher TKE. Overall, the TKE profiles in the porous zone show a distinct difference between the sparse case, where turbulence penetrates vey deep into the bed and the dense case where does not go beyond the centre of the top rod.

Turbulence withi n and Over the Permeable Bed
The turbulent energy production (TKEP) was normalized by dividing the spatially averaged TKEP with ρ.U * 2 ( Figure  15). Fro m the figures it can be concluded that a larger TKEP rate in the dense case compared with the sparse. Further, the larger free stream thickness, have a larger TKEP.
The turbulent dissipation rate (TDR) is the rate at wh ich the turbulent energy is absorbed by breaking the eddies down into smaller until it is ultimately converted into heat by viscous forces. Figure 16 shows the normalized averaged turbulent dissipation rate over the free stream and the porous zone. Fro m the figure, it can be concluded that TDR is lower when free stream thickness is higher, and the porosity effect is smaller than the effect of free stream thickness.

Conclusions
A detailed analysis was made to some flow variables to increase our understanding of turbulent flow through and over permeable bed. These variables were the vertical and longitudinal velocit ies, turbulent kinetic energy, turbulent dissipation rate, and the energy consumed by flow. The study covers free stream, porous and interface zones. Further, the porosity effect and the free stream thickness were also included in the simu lation. The following results were gained fro m the analysis.
• Mean velocities show the expected boundary-layer profile in the free stream zone, and a much slo wer flow in the porous zone. A pipe-like velocity profile fo rms between each two successive rows of cylinders. As expected bed porosity has a major influence on the flow velocities within the porous layer. • TKE in the free stream region have the expected shape with a maximu m value close to the wall and the minimu m close to the free surface. Within the porous layer there is a clear difference between the sparse and the dense arrangement: in the former case TKE penetrates practically throughout the whole porous layer, whereas in the latter case it penetrates only until the half of the top layer of cylinders.
• TKE p roduction and dissipation rate show very similar features. Turbulent shear stress in the free stream shows the expected approximately linear pro file. Within the porous layer, the sparse case has a linear shear stress profile alternating between positive and negative values, typical for a pipe flow, whereas for the dense case the shear stress within the porous layer is practically zero.