1 Introduction

The reduction of turbulent drag in wall-bounded turbulent flows represents one of the great challenges in fluid mechanics, as well as in many industry applications. In these years, the interest in this subject is steadily growing, because even a little bit of reduction can lead to dramatic energy saving in flow systems, such as oil pipelines, turbine blades, and high-speed aircraft wings. Furthermore, the lower fuel consumption results in less pollution and improves environmental sustainability.

The modulation on turbulence by the porous medium is a passive strategy (for which the geometrical equipment is used and no extra energy is inputted) in controlling the flow. The studies in the past have mainly focused on the modulations of the structures and the dynamic of turbulence by isotropic porous surfaces. Zagni and Smith^[1] performed experiments on the open-channel flow over porous beds. They found that the drag and the Reynolds number were higher than those for flows over impermeable walls. They attributed this result to the enhanced energy dissipation caused by the exchange of momentum across the porous interface. Breugem et al.^[2] numerically studied the influence of the highly permeable porous medium on turbulent channel flows. They reported that the turbulent structures near the wall, where the streaks and the associated quasi-streamwise vortices were absent, were quite different from those of the impermeable walls. Moreover, in agreement with other authors, the momentum exchange between the porous wall and the clear fluid was dominated by the larger vortical structures which contributed strongly to the Reynolds stress and drag coefficient. Suga et al.^[3] studied experimentally the effects of wall permeability on the laminar-turbulent transition in channel flows. They observed that the transition to turbulence appeared at lower Reynolds numbers as permeability increased, consistent with the linear analysis by Tilton and Cortelezzi^[4-5]. They also reported that the wall normal fluctuating velocity near the porous wall was enhanced which resulted in the increase in the turbulent shear stress and the higher friction. Rosti et al.^[6] performed numerical studies of turbulent flows over the lower permeable porous wall where the inertial effects were neglected. The porosity and permeability were studied as independent parameters. They illustrated that the permeability played a more important role in modulating the turbulent structures than the porosity. Different from the isotropic porous medium, the studies on the anisotropic porous medium received little attention. Recently, Kuwata and Suga^[7] performed a numerical simulation based on the lattice Boltzmann method. Four kinds of anisotropic porous media were modeled as cubes with square pore arrays aligned along the Cartesian axes, namely, with only vertical permeability, with vertical and streamwise permeabilities, with vertical and spanwise permeabilities, and with vertical, streamwise, and spanwise permeabilities. The simulation was carried out at friction Reynolds numbers of 110 and 230, and revealed that the streamwise permeability had a more significant effect on turbulence than the other directions. Similar to the isotropic porous wall, the turbulent drag was found to be increased in these four anisotropic porous walls. Different from the results in Ref. [7], a drag reduction 20% was obtained by Rosti et al.^[8] in turbulent channel flows over an anisotropic porous wall. They set the permeability in the vertical direction lower than that in the other two directions which limited the wall-normal velocity of the fluid at the interface. In this configuration, the low- and high-speed streak coherence was strongly enhanced which was attributed to the drag reduction. In the present simulation, inspired by the longitudinal riblets, we set the permeability in the spanwise direction lower than that in the other two directions to restrict the spanwise velocity near the wall. We will show that this configuration can also lead to the drag reduction.

In this work, we carry out the direct numerical simulation (DNS) of turbulent open channel flows bounded by a porous wall with the same porosity value 0.6. The flow is driven by the constant pressure gradient which is only applied in the clear fluid domain. Therefore, the friction Reynolds number is fixed at Re_τ=180.

The remainder of the paper is organized as follows. The physical model and mathematical formulation are described in Section 2. For the flow inside the porous medium, we adopt the volume-averaged method^[2] which is governed by the volume-averaged Navier-Stokes equations. The effect of the porous medium is modeled by the Darcy-Brinkman-Forchheimer acting force term^[9-10]. Some statistical quantities of wall-bounded turbulent flows and their validation are given in Section 3. The parametric study and the results are discussed in Section 4. Finally, we summarize our findings in Section 5.

2 Problem and method 2.1 Governing equations and boundary conditions

We consider the turbulent open channel flow as shown in Fig. 1. h is the channel half width. The interface between the porous layer and the fluid is located at y=0, while the impermeable wall is located at y=-h_p. The porosity ε and the permeability K are the parameters that characterize the porous layer. The former is a scalar quantity. The latter is a second-order tensor which reflects the permeability in different directions.

According to the Brinkman-Forchheimer equations^[11], the governing equations are given as follows:

(1)

where ρ is the fluid density, u and p are the volume-averaged velocity and pressure, respectively, ν is the effective kinematic viscosity, and F represents the total force including the resistance force by the porous medium and the driving force^[12-14]. The expression of F is

(2)

where K and B are the permeability and the Forchheimer coefficient tensor, respectively. G is the pressure gradient to maintain the flow. The definitions of K and B are given in Refs. [15] and [16]. In the simulation, we set K_ij=0 if i≠j, K₁₁=K_x, K₂₂=K_y, and K₃₃=K_z. Similarly, B_ij=0 if i≠j. Otherwise, . F_ϵ is the Forchheimer constant. Here, we take a value for the isotropic porous medium given by Ergun^[17] for simplicity,

(3)

where ϵ is the porosity of the porous medium.

It is worth mentioning that, in the fluid domain over a porous medium layer, the porosity can be assumed to be unity, and the permeability is approximately infinite. Therefore, the force F can be seen as G, and Eq. (1) is equivalent to the standard Navier-Stokes equation in this region^[18].

In this study, the periodic boundary conditions are applied in both x- and z-directions with a constant streamwise pressure gradient applied in the clear fluid region. The no-slip and no-penetration conditions are imposed at the bottom wall y=-h_p. The shear-free boundary condition is applied at y=h,

(4)

3 Numerical methods and validation

In order to carry out the present DNS, a fractional-step method developed by Kim and Moin^[19] and Verzicco and Orlandi^[20] is used to solve the governing equations for turbulent open channel flows. Time advancement is performed by the three-stage Runge-Kutta scheme for the convective terms and the Crank-Nicolson scheme with the semi-implicit scheme for the viscous terms. Spatial derivatives are discretized by a second-order central difference scheme.

The computational domain in our simulation is 2πh, h+h_p, πh in the streamwise, normal, and spanwise directions, respectively. The number of grid points is 128× 270× 128 in the streamwise, wall-normal, and spanwise directions, respectively. It is uniform along x- and z-directions with resolutions Δx⁺=8.84 and Δz⁺=4.42, respectively. To resolve the fine turbulent structures near the wall and free surface, a stretched grid which follows a classical coordinate transformation^[21] is applied in the vertical direction with the grid size Δy⁺ ranging from 0.10 to 2.11. It is smaller than the Kolmogorov scale near the wall which means that it can resolve the viscous sublayer and the diffusive sublayer near the boundary.

3.1 Comparison with the DNS of turbulent channel flows

To validate our code, a fully turbulent open channel flow at Re_τ=180 is calculated, and the results are compared with those from computational results and experiments by Handler et al.^[22], Wang et al.^[23], and Komori et al.^[24]. The mean velocity and velocity fluctuations relative to the friction velocity u_τ are shown in Fig. 2. It is seen that our results are in good agreement with the previous results. Although the streamwise velocity fluctuations near the free surface are somewhat lower than those in the references, it can be partly explained that the non-conservative property exists in finite difference schemes as discussed by Morinishi et al.^[25]. The present computational method and relevant code have also been verified in our previous studies^[26-27]. Thus, it can be confirmed that our calculation is reliable for the prediction of the statistical quantities of turbulent open channel flows.

3.2 Comparison with the DNS of turbulent flows with porous walls

Moreover, we carry out the numerical simulation for fully turbulent channel flows with the porous walls at Re_τ=180, and the results are compared with those from the DNS by Rosti et al.^[6]. The computational domain is 2πh, 2.4h, πh in the streamwise, normal, and spanwise directions, respectively. The baseline case with two identical porous walls, whose height h_p=0.2h, is selected as the comparable case. The Darcy number (Da=K/h²) is 1.6×10^-5, and the porosity is ε=0.6. Figure 3(a) compares the mean velocity profiles in the fluid region. The profile is plotted as the difference u₁⁺-U₁⁺ versus the logarithm of the distance from the interface y⁺, where U₁⁺ is the mean velocity of the fluid at the interface. The mean velocity profile by Rosti et al.^[6] is also added as reference, and it turns out to be indistinguishable from the present DNS. We continue the comparison by the root-mean-square velocity fluctuations between the data of Rosti et al.^[6] and the present DNS. As we can see in Fig. 3(b), good agreement is captured in three directions. Thus, it can be confirmed that our code is reliable and precise enough for the prediction of statistical quantities of turbulent channel flows with porous walls.

4 Results and discussion

In our studies, the dimensionless permeability Da=K/h² in the streamwise, vertical, and spanwise directions is set Da_x=Da_y=10^-4 and Da_z=10^-5, respectively. The depth of the porous layer h_p is the only parameter varying from 0.005h to 0.3h. The dimensionless depth h_p⁺=h_p/δ_v∈ (0.9, 54), where δ_v is the viscous lengthscale.

4.1 Mean velocity

Figure 4(a) shows distributions of the mean velocity. The smooth wall case is also added for comparison. We can find that the velocity near the free surface region decreases as the depth of the porous layer increases. However, the lines almost overlap near the interface. Since the flow is driven by a constant pressure gradient for all cases, this indicates that the drag reduction can be obtained by a thinner porous layer, and the drag increase occurs with a thicker porous layer. In the porous region, the velocity quickly decreases to zero regardless of the porous layer thickness, which is attributed to the viscous effect near the wall and the strong resistance in the porous layer.

For the porous wall channel flow, the common way of writing the mean velocity in the log-law region is

(5)

where κ and B⁺ are constant values (κ=0.41, B⁺=5.5 at Re_τ=180). Figure 4(b) shows the mean velocity u⁺-U₁⁺ after removal of the interface velocity versus the distance from the interface y⁺. As we can see, the variation of the velocity profile is negligible in the near-wall region. In the log-law region, the slope of the profile {κ} is unchanged, and the profile seems to substantiate the wall similarity hypothesis by Raupach et al.^[28], which states that for large Reynolds numbers, the turbulence in the outer layer is unaffected by the roughness at the wall. The intercept increases with the decreasing value of h_p⁺ and approaches the smooth wall value B⁺=5.5. This means that the profile monotonically approaches the standard wall law as the effect of porous layer decreases. It is interesting to note that the behavior is consistent with the behavior discovered by Rosti et al.^[6] for the variation of isotropic porous permeability.

4.2 Effects of the porous layer on the flow drag

The drag coefficient is defined as C_f=2τ_w0/(ρU_m²), where τ_w0 is the total shear stress at y=0. As the constant pressure gradient G is applied only in the clear fluid domain in our simulation, τ₀ is constant in different computational cases which equals Gh. The drag reduction ratio^[29] is defined by D_R=(C_f-C_f0)/C_f0, where C_f0 is the wall friction coefficient corresponding to the smooth wall case. Figure 5 shows the measured drag reduction ratio as a function of the dimensionless porous layer depth h_p⁺. As we can see, D_R first decreases at the small porous layer depth, reaches a minimum value and then increases. Two different regimes can be identified, namely, the drag reduction regime and drag increase regime. This is consistent with the behavior we notice in the mean velocity profile. It should be noted that Nikuradse^[30] found that the effect of the rough wall on the drag coefficient was negligible if the height of the roughness element was smaller than the viscous sublayer for the smooth wall case. However, in our present study, D_R < 0 is clearly observed over the range h_p⁺ < 25, and a maximum value (15.3%) is obtained at h_p⁺=9. It is really surprising to note that the optimal depth of the porous layer is consistent with the optimal height of riblets obtained by El-Samni et al.^[31] for turbulent channel flows. However, we obtain a better drag reduction than that obtained with riblets in Refs. [31]-[33] (around 10%). Rosti et al.^[8] also studied the influence of anisotropic permeability on the drag reduction, and they set the vertical permeability Da_y lower than Da_x and Da_z. In their studies, a large h_p⁺ can also lead to the drag reduction. However, in our studies, the thickness of the porous layer is a key parameter that has a significant influence on the turbulent structures and turbulent drag. This discrepancy is related to two main factors, i.e., the mechanism of drag reduction in our simulation is more like riblets which resist the spanwise motion of the fluid near the wall; the vertical permeability in our simulation is much bigger than that in Ref. [8], which allows more turbulent kinetic energies propagating into the porous layer. For the thicker porous layer, the momentum exchange across the interface is enhanced which leads to the unnecessary energy dissipation.

It is known that the drag coefficient in fully developed, incompressible plane channel flows can be divided into a laminar contribution and a turbulent contribution (Fukagata et al.^[34]). When the fluid flows over the porous walls, an extra term occurs which is induced by the porous layer. As derived by Rosti et al.^[6], the total drag coefficient in turbulent channel flows with porous walls is

(6)

where U₁ is the slip velocity at the interface between the clear fluid and the porous layer. U_m is the bulk velocity defined by , and Re=U_mh/ν. Figure 6 shows different contributions to the total drag coefficient as a function of h_p⁺. As we can see, the turbulent contribution C_T increases significantly as h_p⁺ increases, and its distribution shows a parallel trend to that of the drag coefficient C_f. However, the laminar contribution C_L and the porous medium contribution C_P increase a little. This means that the main modification of the drag coefficient by the porous medium is caused by the modulations on the turbulent structures and intensities.

4.3 Turbulence intensities

Results in the coming sections will mainly focus on two distinguished cases of h_p⁺=9 and 36 (the drag reduction case and the drag increase case, respectively). Figure 7 presents the profile of root-mean-square velocity fluctuations which are normalized by the friction velocity at the interface u_τ. The solid line represents the profile of the smooth wall case. As we can see, the significant difference of velocity fluctuations between the drag reduction case and the drag increase case is located near the interface. For the drag reduction case, the velocity fluctuations in three directions decrease compared with the smooth wall case especially for the spanwise velocity fluctuation w_rms. This phenomenon is similar to the velocity fluctuations of channel flows with streamwise riblets^{[31, 35]} and with tensegrity fabrics^[36]. For the drag increase case, v_rms and w_rms increase significantly, and the peak shifts to the interface, but u_rms decreases a lot compared with the smooth wall case. This result is similar to the behavior^{[2, 37]} found in turbulent channel flows with isotropic porous walls. Moreover, the penetration depths of u_rms and v_rms are increased compared with the drag reduction case, while w_rms quickly decreases to zero across the interface. These variations near the interface indicate that the momentum exchange between the porous medium and the clear fluid is enhanced, and the vortical structures lying up the interface have a wider extension especially in vertical and spanwise directions. This behavior is attributed to the weakening of the wall-blocking and wall-induced viscous effects^{[2, 6, 8, 38]}.

As explained in Fig. 6, the contribution of turbulent activity to the drag coefficient is remarkable. According to Eq. (6), the turbulent part is proportional to the averaged shear Reynolds stress. Figure 8 presents the shear Reynolds stress normalized by the bulk velocity U_m. As we can see, the shear Reynolds stress decreases rapidly to zero when moving downwards inside the permeable wall, ~which is a consequence of the resistance of the porous medium. The thinner porous layer reduces the shear Reynolds stress in the clear fluid region, and the thicker porous layer increases it compared with the smooth wall case. As pointed out by Bernard et al.^[39] that the shear Reynolds stress production had a direct relationship to the organized vortical structure. Moreover, the quasi-streamwise vortices were generally observed to lie directly at high skin-friction regions. To demonstrate more vividly, the quasi-streamwise vortical structures distinguished by the Q-criterion^[40-41] are displayed in Fig. 9. It is illustrated that, when the porous layer is thin, the turbulent vortical structures are suppressed which is attributed to the attenuation of the spanwise velocity in the porous layer. As h_p⁺ increases, the number of the quasi-streamwise vortices increases significantly, which indicates that the turbulent intensity is enhanced near the wall. Thus, the momentum exchange between the near-wall region and the outer layer is promoted which sustains a high turbulent drag C_T. This is due to the weakening of the wall-blocking and the wall-induced viscous effects as h_p⁺ increases^{[2, 6, 8, 38]}. It allows the perturbations deep in the porous medium to propagate into the clear fluid region.

To reveal the mechanism of the variations in the quasi-streamwise vortices near the wall, we study the root-mean-square streamwise vorticity fluctuations. Figure 10 displays the streamwise vorticity fluctuations normalized by wall units. The smooth wall case is also included for comparison. It can be seen that the strength of the streamwise vorticity decreases for the thinner porous layer. This behavior is similar to the drag reduction over streamwise riblets, because the porous layer in our simulation mainly impedes the spanwise movement of the fluid in the porous layer. As pointed out by Choi^[42] that the turbulent shear stress in the boundary is produced mainly by the sweeping of the high-momentum fluid, which is a consequence of the legs of hairpin vorticity stretching in the streamwise direction. However, this behavior is reduced as its spanwise movement is restricted. For the thicker porous layer, the strength of the streamwise vorticity is increased, and the center of the streamwise vorticity shifts to the wall. Such behaviors may be attributed to the interaction between the streamwise vorticity living in the clear fluid region and the ejection-like motions arising from the porous layer. This phenomenon is consistent with the bigger velocity fluctuations above the interface for drag increase cases.

Moreover, we perform the quadrant analysis^[43-44] to explore the detailed information about the contributions of different events to the turbulent kinetic energy. Each quadrant is defined as follows. The first Q₁ (u₁' > 0, u₂' > 0) and the third Q₃ (u₁' < 0, u₂' < 0) quadrant events contribute to the positive production of the shear Reynolds stress and the turbulent kinetic energy, respectively. The second Q₂ (u₁' < 0, u₂' > 0) and the fourth Q₄ (u₁' > 0, u₂' < 0) quadrant events correspond to the ejection and sweeping motions, respectively. Near the wall, the major contribution of turbulent kinetic energy comes from the ejection and sweeping events. For the drag reduction case, as shown in Fig. 11(a), the sweeping and ejection events are reduced compared with the smooth wall case, and the other two quadrants change a little. The four quadrants are decreased to zero in the porous layer which are not shown here. The drag increase case is shown in Fig. 11(b). As we can see, the ejection and sweeping events increase significantly compared with the smooth wall case. In the~meantime, the Q₁ and Q₃ events increase as well in the near-wall region. This indicates that there exist quasi-streamwise vortices in the porous layer which can interact with the vortical structures above the interface, which promotes the momentum exchange across the interface.

5 Conclusions

In this work, we mainly focus on the effects of the porous medium with the spanwise permeability smaller than the other two components on the turbulent drag in open channel flows. The flows in the clear region are tracked by the standard DNS, while a volume-averaged approach is applied in the porous region. We show that a drag reduction can be obtained when the dimensionless thickness of the porous layer is less than 25 wall units. The thickness of the porous layer plays an important role in turbulent structures and turbulent intensities. When it is greater than 25, the drag increases compared with the smooth wall case. An optimal drag reduction is 15.3% with the thickness of the porous layer 9 viscous units. An attempt to explore the mechanism for drag variations has been presented by studying different contributions to the drag coefficient, the change in the shear Reynolds stress, the variations of the velocity fluctuations as well as the instantaneous vortical structures. When the porous layer is thinner, it behaves more like a riblet, and the ejection and sweeping events are suppressed due to the restriction of the streamwise vorticity extension in the spanwise direction. When the porous layer is thicker, it behaves more like an isotropic porous medium, and the momentum exchange is promoted between the inner wall and outer wall which enhances the turbulent intensities.

The previous results are obtained only in the numerical simulation. More experiments and DNS with practical design of the solid matrix need to be investigated to confirm the drag reduction effects. In conclusion, our studies may pave a way for the drag reduction with the anisotropic porous wall.