Boundary Layer under Oscillatory Wave

Abstract

Turbulence due to wave motion and propagation is a very important aspect in sediment transport modeling. The boundary layer characteristic during the process will highly influence the sediment transport mechanism at the bottom. 1D model approach has been widely used to assess the turbulent boundary layer. However, the need for a more detailed model leads to the development of a more sophisticated models. This study presents a 2D turbulent model using 𝑘-𝜔 equation to approach the turbulent boundary layer characteristic in oscillatory wave. The equations are solved using finite-difference center scheme. The model is applied to the case of oscillatory sinusoidal and skew wave. The velocity profile and bed stress were compared to those obtained from wind channel experimental data and 1D 𝑘-𝜔 model from previous study.

1. Introduction

The sediment transport phenomenon in coastal area has been studied widely. A detailed and thorough study of this phenomenon requires a better approach in the bed stress approximation. In the near-shore region, the bottom boundary layer characteristic of the wave will play an important role in the sediment transport process. Strong asymmetric shape is often observed in this region due to wave breaking and other coastal process. The boundary layer in this area is usually very thin and, therefore, difficult to observe in nature. For these reasons, experimental and numerical works are preferred to study this phenomenon.

Numerous studies on numerical model have been conducted to understand the turbulent boundary layer under the wave motion. Sana and Tanaka [1] have compared the direct numerical simulation (DNS) data for sinusoidal oscillatory boundary layer over smooth bed to several 1D low Reynolds numbers. A detailed study of the sinusoidal wave profile was given by Suntoyo [2]. He conducted laboratory experiment on oscillatory wave boundary layer using skew and sinusoidal wave profile in wind tunnel. Closed conduits are chosen because they allow generation of high shear stress on the bed. The results are compared to several 1D turbulent models.

1D model was considered adequate for analysis in past researches. However, recent development leads researcher to study in a more detailed process of wave propagation and run-up. Convective acceleration plays an important role during the wave propagation and run up in shallow water. Therefore, the 1D model is no longer adequate to assess the boundary layer. The 1D model neglects the convective acceleration effects in the flow. The importance of the convective acceleration has been shown by Adityawan et al. [3]. They had made comparison between 2D model with convective acceleration and no convective acceleration for sinusoidal wave. The results showed that the convective acceleration term should not be neglected.

There are several types of models for turbulence flow. The two equation models are commonly used due to its reliability. The 𝑘-𝜀 model was one of the first two equation models used to study turbulence. 𝑘-𝜔 model is a further development from the 𝑘-𝜀 model. In this study, a 2D vertical model is developed based on the Reynolds-averaged equation.

The present study will focus on the investigation of turbulent boundary layer characteristic using 2D model. The model is developed using the Reynolds-averaged equation and 𝑘-𝜔 turbulence model. The model is applied to the experimental studies and 1D numerical model of oscillatory wave from previous study. This study will focus on the case of oscillatory sinusoidal wave and skew wave. The velocity profile and the bottom shear stress from 2D calculation are compared to the experimental data and 1D calculation from previous study.

2. 𝑘-𝜔 Model

The governing equation for the model is based on the Reynolds-averaged equations of continuity and momentum as follows: 𝜕𝑈_𝑖𝜕𝑥=0,(2.1)𝜌𝜕𝑈_𝑖𝜕𝑡+𝜌𝑈_𝑗𝜕𝑈_𝑖𝜕𝑥_𝑗=−𝜕𝑃𝜕𝑥_𝑖+𝜕𝜕𝑥_𝑖2𝜇𝑆_𝑖𝑗−𝜌𝑢^_𝑖𝑢^_𝑗,(2.2) where 𝑈_𝑖 and 𝑥_𝑖 denote the mean velocity and location in the grid, 𝑢^_𝑖 is the fluctuating velocity in the 𝑥(𝑖=1) and 𝑧(𝑖=2) directions, 𝑃 is pressure, with the pressure gradient itself being a function of the free stream acceleration, 𝜈 is the kinematics viscosity, 𝜌 is the density of the fluid, 𝑆_𝑖𝑗 is the strain-rate tensor, and −𝜌𝑢^_𝑖𝑢^_𝑗 is the Reynolds stress tensor.

We have the following: 𝑆_𝑖𝑗=12𝜕𝑈_𝑖𝜕𝑥_𝑗+𝜕𝑈_𝑗𝜕𝑥_𝑖,(2.3) where −𝜌𝑢^_𝑖𝑢^_𝑗 is the Reynolds stress tensor. The Reynolds stress tensor is given through eddy viscosity by the Boussinesq approximation:−𝑢^_𝑖𝑢^_𝑗=𝑣_𝑡𝜕𝑈_𝑖𝜕𝑥_𝑗+𝜕𝑈_𝑗𝜕𝑥_𝑖−23𝑘𝛿_𝑖𝑗(2.4) with 𝑘 being the turbulent kinetic energy and 𝛿_𝑖𝑗 being the Kronecker delta.

The turbulence model is the 𝑘-𝜔 model. 𝜔 is the ratio of turbulence dissipation rate 𝜀 to the turbulent kinetic energy 𝑘 and may be regarded as a characteristic time scale of the turbulence. The 𝑘-𝜔 model equations are with the following equation: 𝜕𝑘𝜕𝑡+𝑈_𝑗𝜕𝑘𝜕𝑥_𝑗=𝜏_𝑖𝑗𝜕𝑈_𝑖𝜕𝑥_𝑗−𝛽^∗𝑘𝜔+𝜕𝜕𝑥_𝑗𝑣+𝜎^∗𝑣_𝑡𝜕𝑘𝜕𝑥_𝑗,𝜕𝜔𝜕𝑡+𝑈_𝑗𝜕𝜔𝜕𝑥_𝑗=𝛼𝜔𝑘𝜏_𝑖𝑗𝜕𝑈_𝑖𝜕𝑥_𝑗−𝛽𝜔²+𝜕𝜕𝑥_𝑗𝑣+𝜎𝑣_𝑡𝜕𝜔𝜕𝑥_𝑗,𝑣_𝑡=𝑘𝜔.(2.5) The values of the closure coefficients are given by Wilcox [4]. They are𝛽=3/40,𝛽^∗=0.09, 𝛼=5/9, and 𝜎=𝜎^∗=0.5.

The use of 2D model applies that convective acceleration is not neglected in the penetrating pressure gradient, coming from beyond the boundary layer. The pressure gradient at the upper boundary layer is calculated as follows: 𝜌𝜕𝑈_𝑖𝜕𝑡+𝜌𝑈_𝑗𝜕𝑈_𝑖𝜕𝑥_𝑗=−𝜕𝑃𝜕𝑥_𝑖.(2.6) Given that the gradient in vertical direction in the free stream boundary is equal to zero, then (2.6) becomes 𝜌𝜕𝑈_𝑖𝜕𝑡+𝜌𝑈_𝑗𝜕𝑈_𝑖𝜕𝑥_𝑗=−𝜕𝑃𝜕𝑥_𝑖(2.7) in which 𝑖=1. No slip boundary condition is applied at the bottom layer; hence 𝑢,𝑣,𝑘=0. The advantage of k-𝜔 model as compared to the 𝑘-𝜀 is the introduction of wall function to impose the value of 𝜔 directly at the wall. Bed boundary condition is given for 𝜔 value at the wall directly, given as follows: 𝜔_𝑤=𝑈^∗𝑆_𝑅/𝜈,(2.8) where 𝜔_𝑤 is the specific dissipation rate at bottom, 𝑈^∗ is friction velocity or (𝜏₀/𝜌) and parameter 𝑆_𝑅 is related to grain roughness Reynolds number 𝑘⁺_𝑠=𝑘_𝑠𝑈^∗/𝑣, with 𝑘_𝑠 corresponding to roughness. 𝑆_𝑅 can be calculated as follows: 𝑆_𝑅=50𝑘⁺_𝑠²,𝑘⁺_𝑠<25,𝑆_𝑅=100𝑘⁺_𝑠,𝑘⁺_𝑠≥25.(2.9) Smooth wall can be approached as 𝑘⁺_𝑠<5 for stability. The governing equations above are solved using finite-difference scheme, in which the horizontal axis component is solved directly using finite-difference center scheme where as the vertical axis component is solved using finite-difference Crank Nicholson scheme on a staggered grid. For example, rearranging (2.2) in 𝑥 direction will provide the following: 𝜕𝑢𝜕𝑡+𝑢𝜕𝑢𝜕𝑥+𝑣𝜕𝑢𝜕𝑧=−𝜕𝑃𝜕𝑥+𝜕𝜕𝑥2𝜈𝑆_𝑥𝑥−𝑢^_𝑥𝑢^_𝑥+𝜕𝜕𝑧2𝜈𝑆_𝑥𝑧−𝑢^_𝑥𝑢^_𝑧,𝑢𝜕𝑢𝜕𝑥−𝜕𝜕𝑥2𝜇𝑆_𝑥𝑥−𝑢^_𝑥𝑢^_𝑥−𝜈𝜕²𝑣𝜕𝑥𝜕𝑦−𝜈_𝑡𝜕²𝑣𝜕𝑥𝜕𝑦+𝑣𝜕𝑢𝜕𝑦_{spatial-solvedusingfinite-differenceforwarddifferencescheme}+𝜕𝑃𝜕𝑥_{backwarddifferencefromfreestreamtreatedasconstantintheCrankNicholsonmatrix}+𝜕𝑢𝜕𝑡=𝜈𝜕²𝑢𝜕𝑦²+𝜈_𝑡𝜕²𝑢𝜕𝑦²_{spatial-solvedusingfinite-differencestaggered.CrankNicholson}(2.10)

3. Simulation and Results

3.1. Case Study

The model is used to investigate the boundary layer characteristic under oscillatory sinusoidal and skew wave experiment [2]. The upper boundary condition will be given as the pressure gradient according to the sinusoidal velocity profile of the flow. The wave parameters for the sinusoidal and skew wave profiles are as given in the experiment as shown in Figure 1 and Table 1.

Table 1 

Wave parameters [2].

Figure 1 

Wave profile [2].

The 2D model requires a new definition of the model domain. Therefore, the model domain is adjusted to be equal to the wave length. The wave length is approached by shallow water wave theory. This condition is chosen because the experiment was conducted under the assumption of a very long wave.

Periodic boundary condition is applied to the left and right free end of the domain (Figure 2). The boundary is applied so that the parameters at grid 1 would be equal to the parameters at grid 𝑖𝑚-1 and the parameters at 𝑖𝑚 would be equal to the parameters at grid 1, with 𝑖𝑚 corresponding to maximum number of grids. No slip boundary condition is applied at the bottom layer.

Figure 2 

Periodic boundary.

The grid spacing for 𝑥 axes is 0.05 𝑧_ℎ, and the grid spacing for 𝑦 axis is starting at 0.005 𝑧_ℎ with grid refinement to the free stream, in total number of 100. Here, 𝑧_ℎ is the distance from wall to the free stream. The smaller grid spacing in 𝑧 axis is required to capture the boundary layer profile.

3.2. Velocity Profile

Comparison was made to the velocity profile and bottom shear stress obtained from the 2D model to the parameters obtained from experiment and 1D model for both scenarios. Observation grids are located in the center line of the model domain. The vertical flow distribution is shown for every 1/8 of the wave period (𝑇/8). The bottom shear stress is shown as a function of time.

The simulation result for velocity profile from the simulation of skew wave is shown in Figure 3(a). The model is able to simulate the turbulent boundary layer due to skew wave. The result is compared to the velocity profile obtained from experiment and 1D 𝑘-𝜔 by Suntoyo [2]. The present model shows good comparison with experimental data. It also shows a very good match to the 1D model result.

(a)

(b)

Figure 3 

Velocity profile results comparison; experiment data [2], 1D 𝑘-𝜔 [2], 2D 𝑘-𝜔; skew wave (a), sinusoidal wave (b).

The simulation result for velocity profile from the simulation of sinusoidal wave is shown in Figure 3(b). The model is able to simulate the turbulent boundary layer due to sinusoidal wave. The result is compared to the velocity profile obtained from experiment and 1D 𝑘-𝜔 by Suntoyo [2]. As in the skew wave scenario, this scenario also showed good comparison to the experiment data and the 1D model. It is also observed that the 2D model tends to overestimate the velocity during the rundown. In general, the model results show a good comparison to the 1D model. The 2D model acts similarly to the 1D model.

3.3. Bed Shear Stress

The second parameter observed is the bottom shear stress. The bottom shear stress obtained from the model is compared to the empirical formulas, 1D model, and experimental data. Several empirical formulas are also used for comparison.

The bottom shear stress in the experiment [2] was estimated by fitting the logarithmic velocity distribution to the measured velocity𝑢=𝑈^∗𝜅ln𝑧𝑧₀,(3.1)𝜏₀=𝜌𝑈^∗(𝑡)||𝑈^∗(𝑡)||,(3.2) where 𝜅 is the von Karman’s coefficient (0.4) and 𝑧₀ is the value of 𝑧 at which the logarithmic velocity profile predicts zero velocity.

In this study, two empirical methods for estimating bed stress under the skew wave profile are used for comparison. Kabiling and Sato [5] estimated the bed shear stress within a basic harmonic wave cycle modified by the phase difference. The bed stress is proportional to the square of stream velocity𝜏₀𝑡−𝜑2=12𝜌𝑓_𝑤𝑈(𝑡)||𝑈(𝑡)||.(3.3) Here, 𝜑 is the phase difference and 𝑓_𝑤 is the wave friction factor.

Nielsen [6] suggested that the bed stress is proportional to the square of instantaneous friction velocity incorporating the acceleration effect under asymmetric wave𝑈^∗(𝑡)=𝑓_𝑤2𝜑𝑈(𝑡)+sin𝜑𝜎𝜕𝑈(𝑡)𝜕𝑡.(3.4) The bed shear stress is obtained by applying (3.4) to (3.2).

The time variation of wave-induced bottom shear stress is characterized by a large peak over a very short time interval preceding the wave crest. Both 1D and 2D turbulent boundary layer models are superior to empirical formula. The bed stress profile from these models shows better comparison to the experiment data compared to the empirical formula. The comparison for the scenario of skew wave is given in Figure 4(a) and for the scenario of sinusoidal wave is given in Figure 4(b).

(a)

(b)

Figure 4 

Bed shear stress comparison with 2D 𝑘-𝜔; skew wave (a), sinusoidal wave (b).

Both scenarios give good comparison to shear stress from calculation and experiment. The model tends to give overestimation of the shear stress magnitude, higher than the 1D model or experiment in closed conduits. The effect of convective acceleration is neglected in both 1D model and closed conduit experiments.

4. Conclusion

The 2D 𝑘-𝜔 model has been used to investigate the boundary layer characteristic under the oscillatory wave profile over smooth surface. Two scenarios are chosen, skew wave profile scenario, and sinusoidal wave profile scenario. Two parameters were observed in the simulation: the velocity profile and the bed shear stress. Velocity and bed shear stress in the boundary layer obtained from the model has shown good comparison with both experimental and calculation result. Overall, the 2D model results for both scenarios show similar performance to the 1D model. However, it tends to provide a higher bed stress. The 2D model assessed the pressure gradient term as a function of local acceleration and convective acceleration. The convective acceleration effect is not accounted in the 1D model or oscillatory experiment.

The proposed model has been implemented to assess boundary layer properties on flat bed showing that convective acceleration may have important role in the bed stress production. Further investigation should be conducted with respect to open free surface on sloping bed.

Acknowledgments

This research was partially supported by Grant-in-Aid for Scientific Research from the Japan Society for Promotion of Science (no. 21360230, no. 22360193), the River Environmental Fund (REF) from the Foundation of River and Watershed Environmental Management (FOREM), Japan, and Open Fund for Scientific Research from the State Key Laboratory of Hydraulics and Mountain River Engineering, Sichuan University, China. The first author is a scholarship holder under the Indonesian Ministry of Education auspice.