A NUMERICAL STUDY OF WAVE-CURRENT INTERACTION IN THE BOTTOM BOUNDARY LAYER

In the present work, a numerical wave-current flume has been developed, based on a standard k-ε model. The numerical flume was 12.86m in length, with a numerical beach at one end of the flume. The Volume of Fluid (VOF) method was used to capture the free surface in the flume. The velocity profile obtained at the test section from the numerical simulation has then been compared with experimental data and good agreement found. Periodic velocities in the bottom boundary layer have been obtained which agree well with the experimental data. The model provides an insight to the changes in bed shear stress time histories that characterise wave current interaction.


INTRODUCTION
The hydrodynamics of coastal zones are dominated by two important factors: gravity waves and currents.The waves are usually generated by the wind, while the currents can be generated by the waves, density variations, tides etc.With the development in marine energy exploitation, the interaction between surface gravity waves and a tidal current (WCI) is highly significant for engineers concerned with hydrodynamics in the near-shore region.Renewable energy devices such as tidal turbines, as well as other forms of coastal structure, are located in a complex environment involving ocean waves and a tidal current.Evaluation of coastal erosion, design of harbour structures, pipelines and tidal turbines, and dispersal of pollutants are examples which require a clear understanding of wave-current interaction.Physically, it is complicated to determine how ocean waves interact with the turbulent currents because of many complex physical processes involving various temporal and spatial scales.Considerable research has been carried out to understand the dynamics of wave-current interaction in the past 60 years.An enhanced bed shear stress together with altered velocity distributions are important phenomena in a combined condition.An eddy viscosity concept has been adopted widely by many researchers, and various theories on eddy viscosity distributions with space and/or time have been developed.This was first presented by Lundgren (1972), based on the vectorial combination of a current eddy viscosity and wave eddy viscosity.This is only applicable to a wave-dominated condition because the characteristics of the wave boundary layer were determined from the wave-alone case, without considering the current effects.Grant and Madsen (1979) parameterized the effect of waves on a turbulent current by an enhanced apparent bottom roughness near the bed, based on a time-invariant turbulence eddy viscosity model.This was later modified by Malarkey and Davies (1998), to include time variations in the eddy viscosity.Christoffersen and Jonsson (1985) developed different models applicable to two roughnesses, using different time-independent eddy viscosity distributions.Myrhaug and Slaattelid (1990) investigated bottom friction coefficients and the phase lead of the bottom shear stress over the free-stream oscillatory velocity, relying on a similarity theory and a time-invariant eddy viscosity distribution.A mixing length turbulence model was applied by Bakker and von Doorn (1978) to investigate the bottom boundary layer, taking into account temporal variations inside the wave boundary layer.Fredsoe (1984) used the momentum equation to study the bottom boundary layer.It was assumed that time scales for the change in the outer velocity are much larger than that for the decay of eddies formed in the wave boundary layer.From this they deducted that the theory breaks down at a very high wave frequency, where the history effect of turbulence formed in the previous half-cycle is no longer negligible.A review of these models is included in a review paper by Deigaard et al. (1985), who further evaluated the corresponding application to pipeline stability.Soulsby et al. (1993) reviewed several analytical models and gave a parameterized model by curve fitting.But time series of bed shear stress cannot be obtained from this model, and therefore the turbulent memory effects of reversing flow during previous half-cycles cannot be obtained from these models.Previous experiments on WCI in a vertical plane (Kemp andSimons, 1982, 1983;Klopman, 1994;Umeyama, 2005Umeyama, , 2009) ) have shown that the logarithmic velocity profile, which is characteristic of a unidirectional turbulent flow, is changed with the addition of wave motion.This is particularly obvious in the vicinity of the bed and free surface.Many theories have been proposed to explain this phenomenon, mainly based on the waveinduced stresses (Nielsen and You, 1996;You, 1996;Groeneweg and Klopman, 1998) or secondary circulations in the cross section (Nepf and Monismith, 1991;Dingemans et al., 1996;Klopman, 1997;Groeneweg and Battjes, 2003).
Past studies have identified the non-linear characteristics of WCI, with experimental studies describing the phenomena and semi-empirical models attempting to quantify them.This paper describes a CFD model to provide a better understanding of what is happening in the bottom boundary layer under a combined wave and current condition.Numerical results of boundary layer dynamics over a smooth bed are presented.

METHODOLOGY
A numerical wave-current flume has been developed using a CFD package (Ansys Fluent), based on the finite-volume scheme.The Reynolds-averaged Navier-Stokes equations were adopted for modelling the incompressible flow and the Volume of Fluid method was used to track the free surface.A steady turbulent current boundary layer was developed from a prescribed uniform velocity at the inlet, using a standard k-ε model.In order to fully resolve the bottom boundary layer, a Low-Reynoldsnumber scheme was adopted instead of the High-Reynolds-number scheme adopted before (Teles et al., 2013).Second-order Stokes waves were generated through the velocity inlet boundary condition.An artificial damping zone was placed near the pressure outlet boundary to avoid wave reflection.The experimental data from Kemp and Simons (1982) were used in order to verify the numerical results.

Governing Equations
The numerical wave-current flume was based on the RANS equations and an enhanced near-wall treatment was adopted to fully resolve flow field in the vicinity of the bed (Ansys, 2013).This is a nearwall modelling method specifically appropriate for a very refined mesh near the bed (typically with the first near-wall node placed at y + ≈1).

Initial and Boundary Conditions
In order to generate an initial turbulent current, similar boundary conditions to those for the case 'CA' (Kemp and Simons, 1982) were adopted here.This involves a pressure-inlet boundary condition, with streamwise velocity U of 0.2m/s and turbulence parameters (turbulence intensity of 3.9% and hydraulic diameter of 0.43m, both the same as in the 'CA' case).An inlet streamwise velocity U of 0.2m/s was chosen, larger than the experimental one (0.185m/s) since a pressure-inlet boundary condition was specified at the inlet instead of a velocity-inlet boundary condition.Different velocity profiles have also been tried during preliminary tests by the authors (by user-defined functions of velocity profiles), in order to find the most appropriate boundary conditions for the model.The outlet boundary condition was a pressure-outlet, specifying the same turbulence parameters as the inlet specification.Free surface level (0.2m) and bottom level (0m) were both specified through the pressureinlet and pressure-outlet boundary condition, which is important in a two-phase multiphase flow problem.The atmospheric pressure was specified in the model as 1atm = 101325 Pa in the air domain.The smooth bed was again specified as a non-slip wall, without any roughness involved.The top boundary of the computational domain was specified as a slip-wall, with zero shear stress.After a stable current was obtained, waves were added onto the current.In this case, the wave condition 'WA1' was superimposed onto the background current after a flow time of 19.6s.The characteristics of 'WA1' were input based on the Doppler Shift theory (Peregrine, 1976) to ensure that the observed frequency is unchanged (T =1s) under a combined condition.After the initial turbulent current was generated, some changes in boundary conditions and cell zones were made to superimpose wave motions.The waves were generated through the velocity-inlet boundary condition, specifying the current velocity magnitude of 0.185 m/s and turbulence parameters (turbulent intensity and hydraulic diameter remained the same as 'CA').Several other current velocity magnitudes were also tried by the authors: the magnitude of 0.185 m/s gave the closest results compared with experimental data.Wave parameters were specified at the inlet, through a wave boundary condition.Second-order Stokes waves were generated, with a wave height of 0.02m, wave length of 1.446m, phase difference and wave heading angle of 0 as used in the experimental tests.A numerical beach was then switched on, covering the region from x=8.52m to the end of the flume in order to avoid wave reflections.This is accomplished by adding a damping sink term to the right-hand of the momentum equation for cells in the zone near the pressure-outlet boundary.The damping resistance was set to 100 (m -1 ) in order to suppress the reflected waves.

Description of solve-settings
A transient type pressure-based solver was adopted in the generation of a steady current.In terms of spatial discretisation, second-order accuracy is necessary and therefore adopted in the model for its high resolution.The implicit volume fraction scheme was adopted to solve the phase continuity equation iteratively together with the momentum and pressure.Residual target was set as 0.00001 for the convergence criteria, which was achieved within each time step.The time stepping method was chosen as fixed, with 0.001s chosen as the time step size for its greater accuracy and computational speed.For the pressure interpolation scheme, PREssure STaggering Option (PRESTO!) algorithm has been used as recommended for all VOF calculations (Ansys, 2013).This scheme uses the discrete continuity balance for a 'staggered' control volume about the face to compute the face pressure.The interface scheme 'compressive' is adopted as an interpolation scheme for its high accuracy and medium-to-high computational speed (Ansys, 2013).This is a second-order reconstruction scheme.In terms of temporal discretisation, a second-order implicit transient formulation is adopted in the model for its high accuracy.

RESULTS AND DISCUSSIONS
The mean horizontal velocity profiles under a combined condition, on linear and semi-logarithmic scales, are shown in Figure 2(a) and 2(b) respectively.Good agreement is found with the laboratory data from Kemp and Simons (1982), with the maximum relative error within 3%.The periodic velocities near the smooth bed, at 18°phase intervals through half a wave cycle, are shown in Figure 3.This also shows good agreement with the experimental data.A convergence study of the time step size and mesh has been conducted in the present study, which has not been included in any previous work.The basic case adopted above (50000 cells and time step size of 0.001s) was chosen to present results after the convergence tests were performed.This demonstrated that the solution converged further with a decreasing time step size and increasing number of mesh cells (Figure 4 and 5).Due to the high computational cost, the number of mesh cells was limited to 200000 here and no more refined mesh is discussed.Figure 6 gives the velocity distributions at 18°intervals within full depth.It is found that flow reversal takes place very near the bed, within about 0.4mm in the combined wave-current condition.This is in good agreement with the experimental data of Kemp and Simons (1982).The flow reversal involves an inflexion point, which is normally a sign of flow instability.
Figure 7(a) and (b) show the time series of surface elevation and bed shear stress in a combined wave-current condition.This provides more insights regarding the bed shear stress time histories, which is important for sediment transport.The figure demonstrates the phase advance between the free surface elevation and the bed shear stress, confirming that the peak of bed shear stress occurs prior to the wave crest.

SUMMARY
A CFD model has been developed using standard software to simulate flows with combined waves and currents.The turbulent current was generated by a pressure-inlet boundary condition and a pressure-outlet condition at the other end.The expected wave was generated by giving a velocity-inlet boundary condition.Reflected waves were absorbed by the numerical beach, which is based on adding a sink term into the momentum equation.A convergence test was carried out and the basic test (50k mesh cells, time step size of 0.001s) takes 60 hours, which is practical for numerical modelling.Results show good agreement with experimental data over a smooth bed for mean velocities and periodic velocities.Flow reversal was observed in the vicinity of the bed.
The results of bed shear stress time histories indicate a phase lead compared with the free surface elevation.The results can be useful for wind farms designers and sediment transport studies, where the time series of bed shear stress are highly significant.It also suggests further research to investigate equivalent phenomenon under real wave environments, either by numerical experiments or by field tests.

Figure 1 .
Figure 1.-Model grid.Part of the mesh, combined wave-current flow

Figure 3 .
Figure 3. -Periodic velocities in the bottom boundary layer, WCA1 (CFD results presented as lines, experimental data presented in dots).

Figure 4 .
Figure 4. -Convergence of numerical results for mean horizontal velocity profile for the case WCA1.Meshes with 50k cells, and two time steps = 0.005s and 0.001s.

Figure 5 .
Figure 5. -Convergence of numerical results for mean horizontal velocity profile for the case WCA1.Time step = 0.001s, mesh size = 50k and 200k cells.