THE INFLUENCE OF THE TURBULENCE CLOSURE MODEL ON WAVE-CURRENT INTERACTION MODELING AT A LOCAL SCALE

An advanced CFD solver based on the RANS (Reynolds Averaged Navier-Stokes) equations is used to evaluate wave-current interactions through numerical simulations of combined wave-current free surface turbulent flows. The repercussions of various schemes for modeling turbulence effects is addressed with a special attention to the exchanges and fluxes of momentum and energy between the mean flow components and the wave (oscillatory) component. Numerical simulations are compared with experimental data from Klopman (1994).


INTRODUCTION
The co-existence of a variety of phenomena, with a wide range of temporal and spatial scales, makes the study and the modeling of the coastal zone a difficult task.The tides and the intrinsic associated modulation, the high variability of the wind forcing, the irregularity of the bathymetry and of the littoral morphologies, the turbulent fluctuations and the inhomogeneity of the different current fields have all to be taken in to consideration.On the other hand, the knowledge of the existing physical processes is essential, namely, for the design of maritime structures, the optimization of renewable energy potential or the understanding of the littoral morphodynamics.
The study of the combined effect of waves and currents in free surface flows is of relevance and a major issue in nearshore hydrodynamics.Aiming to shed light on the subject and increase knowledge of the fundamentals of this interaction, a number of laboratory experiences had been designed.Among the most important, lies Klopman (1994), who carried out laboratory experiments in a wave flume with a rough bed.He obtained vertical profiles of mean horizontal velocity and amplitude for three distinct cases: only waves, waves following the current, and waves opposing the current.A test with only a mean constant discharge was also prepared, from which it was measured mean horizontal velocity and shear stress' vertical profiles.Klopman observed that when waves were following the current the mean horizontal velocity increased from the bed until a certain level and then started to decrease.On the other hand, when waves were opposing the current the mean horizontal velocity increased until the free surface.Similar results were obtained by Kemp andSimons (1982, 1983) and Umeyama (2005).
At the same time, numerical solutions have been developed to evaluate the combined effects of waves and current.Dingemans et al. (1996) proposed the inclusion of the Craik-Leibovich vortex force in the mean-current equations.Groeneweg and Klopman (1998) and Groeneweg and Battjes (2003) presented a Generalized Lagrangian Mean approach (GLM), and Olabarrieta et al. (2010) used a three dimensional Navier-Stokes equations model.
The main objective of the present work is to study the interactions between the surface waves and the currents, with emphasis on the changes of the mean horizontal velocity profile induced by the wave's propagation on the free surface.
To accomplish this purpose, it was used the CFD (Computational Fluid Dynamics) model Code_Saturne (Archambeau et al. 2004) that solves the primitive equations of Newtonian flows with the Reynolds decomposition for the mean flow and turbulent fluctuations.
The model was developed at EDF R&D, Chatou, France and was initially designed for pressurized flows in industrial facilities and networks.Subsequently, the model was adapted for the study of the interaction between waves and currents considering turbulence effects on free surface flows.For dealing with the free surface, the Arbitrary Lagrangian Eulerian (ALE) methodology (Archambeau et al. 1999) was used.
The Reynolds Averaged Navier-Stokes (RANS) equations need a hypothesis to close the system of equations.Among the turbulence closure models available in the code, it was made the choice of using the first order two-equation models, k−ε and k−ω, largely applied for their simplicity, and the second order Reynolds stress transport model R ij −ε.The choice of the turbulence closure model has a major

Introduction
The numerical model Code_Saturne (Archambeau et al. 2004), a RANS solver was used for the simulation of waves and current interactions.Code_Saturne was developed at Electricité de France (EDF) and it is a general purpose CFD open source software (available at http://www.codesaturne.org).It solves the Reynolds-averaged Navier Stokes equations for the mean incompressible flow by employing a finite-volume discretization approach.
The equations are written in a conservative form and then integrated over control volumes.
( ) 0 div u t S u represents additional momentum source terms that can be prescribed by the user, ρ the fluid density, u the mean flow fluid velocity vector and σ the stress tensor.For a turbulent flow, the stress tensor σ includes the effect of pressure, viscous stress τ, and the turbulent Reynolds stress tensor R ij .
To represent the free surface variation, Code_Saturne uses a moving-mesh approach, the Arbitrary Lagrangian Eulerian (ALE) methodology (Archambeau et al. 1999).With this module, the Navier Stokes equations gain a new term which involves the vertical velocity of the mesh.For each time step, the mesh is updated accordingly.

Turbulence modelling
In Code_Saturne it is possible to choose a number of first and second order turbulence closure models to model the Reynolds stress tensor R ij and close the system of equations ( 1) and (2).While in a first order turbulence model the Reynolds stress tensor is linked to the mean flow velocity through the Boussinesq hypothesis and the turbulent viscosity approximation, in a second order turbulence model the Reynolds stresses are solved explicitly with transport type equations.
In the present work, it was chosen to use the first order two equation models k−ε with linear production (Guimet and Laurence 2002) and the k−ω SST (Menter 1994), widely used for their simplicity, and the second order Reynolds stress transport model R ij −ε SSG (Speziale et al. 1991).
k−ε with linear production (Guimet and Laurence 2002) and R ij −ε SSG (Speziale et al. 1991) are the so called High Reynolds models, for which one should guarantee that the thickness of the first computational cell near the wall is larger than the thickness of the viscous sublayer.Consequently, an analytical treatment (wall functions) is needed in the near wall area.On the contrary, the k−ω SST (Shear Stress Transport), proposed by Menter (1994), exhibits a better behavior near the wall.Therefore, no wall functions are needed but also the thinner mesh is confined to the near wall zone.
In the first order turbulence models additional equations are needed to calculate the turbulent viscosity.The transport of the turbulent kinetic energy k (3) and turbulent dissipation rate ε (4) are examples of equations commonly used as closure, leading to the well-known k−ε turbulence model. .
k is the turbulent kinetic energy, ε the turbulent dissipation, µ the fluid's dynamic molecular viscosity, µ t the turbulent viscosity, P accounts for the production of the kinetic energy through mean shear stresses, G the production term related to gravity effects and finally In the k−ω SST model (Menter 1994), equation ( 3) is solved for k, but for the dissipation a so called specific turbulence dissipation ω = ε/k is used.
The Reynolds Stress Models (RSM) do not include the eddy viscosity hypothesis and the Reynolds stress transport equation ( 6) accounts for the directional effects of the Reynolds stress fields and.In RSM there are six transport equations for the six components of the Reynolds stress tensor and one equation for the dissipation rate (7).
. ( ) R ij represents the Reynolds stress tensor, P ij and G ij are the turbulence production tensors related to mean shear stresses and gravity effects respectively, Φ ij is the pressure strain term, d ij and d ε the turbulent diffusion terms and ε ij the dissipation term (considered isotropic) (Archambeau et al. 2004).

Model setup
The mesh generation is subjected to a number of conditions that the modeller has to take into account.
As referred, High Reynolds models were used and some constraints appear relatively to the size of the neighbourhood cells of the bottom.At the same time, and because some important effects were meant to be analyzed in this region, such as the roughness' influence on the vertical profile of the measured quantities, a high resolution was required.
To guarantee a good representation of wave's propagation along the channel, the mesh was discretized in about 10 cells per wave length.Right next to the moving wall (wave maker) the mesh could not be too refined to avoid mesh crossover and cause the divergence of the simulation.The representation of the generated mesh is given on Figure 1.
To generate the numerical waves, a second order sign, which minimize undesirable free superharmonics and sub-harmonic waves, was imposed on the lateral moving wall.Waves propagated in positive x-direction.The following expression for the wave board motion X 0 (t) in equation ( 8) (Dean and Dalrymple 1991) was introduced: K represents the wave number, h the water depth, H the wave height, σ the relative angular wave frequency and t the time.To avoid a sudden mesh's horizontal movement and thus mesh crossover, the signal given by (8) had to be progressively imposed at the lateral boundary.To minimize the wave's reflections at the end of the numerical channel, the latter was extended and filled with an artificial fluid that has a viscosity which is progressively increasing (Figure 2).Recently, Cozzi (2010) has adapted the ALE method for representing wave propagation on free surface flows for a ideal fluid.One of the conditions imposed at the free surface, to guarantee a zero net mass flux, is represented by the following expression, . .w represents the vertical velocity of the mesh referred above, m fS is the mass flux in the case of considering the free surface with a fixed mesh, S z the vertical component of the normal vector to the free surface.
In the present study, the purpose was to study wave current interactions with Code_Saturne taking into account turbulence effects on free surface flows.It was verified that an additional condition (11), proposed by Celik and Rodi (1984), had to be imposed on the free surface.This condition was essential to get the right turbulence effects, since it accounts for the reduction of the length scale of turbulence near the free surface.As discussed by Nezu and Nakagawa (1993), with this boundary condition, the turbulent dissipation, which determines the length scale, will be higher than the value that would follow from a zero-gradient condition.Hence from equation ( 5) it can be seen that the eddy viscosity decreases towards the free surface.ε and k are the values of the turbulent dissipation and turbulent kinetic energy at the water surface, respectively, α = 0.18 is an empirical constant and h is the water depth.The RSM model does not compute explicitly the turbulent energy k, it is estimated as half the sum of the normal Reynolds stresses (12).

(
) LABORATORY DATA To validate the numerical results, experimental data obtained by Klopman (1994) was used.Klopman (1994) carried out a series of laboratory experiments in a channel roughened with sand with 46 m long (x direction), 1 m wide (y direction) and water depth established at h = 0.5 m.The channel was provided with a recirculation system with a constant discharge of about Q ≈ 80 l/s.
The channel had two computer controlled wave boards, one to generate the waves with a second order signal to minimize free long waves, and another one to absorb the waves.
Numerical simulations were carried out for the test cases with current and waves only, and monochromatic waves following and opposing current.The wave height was H = 0.12 m and the period T = 1.44 s.Mean horizontal velocity profiles and horizontal velocity amplitudes were measured by a laser-Doppler velocimeter for each test at the middle of the channel (x = 22.5 m and y = 0.5 m).For the case of only current a description of the shear stress was also made.

Case with only a current
It should be emphasized that this sensitivity test is made with the default parameterization set of the Code_Saturne turbulence modeling.This means that no optimization was tried with the free parameters of each turbulence closure model.Also the boundary condition on free surface (11) was imposed for each model.
The first case to be considered in Code_Saturne was the imposition of a mean discharge in the channel, without generating waves.In Figure 3, the mean horizontal velocity profiles from numerical results using the three turbulence closure models: k−ε, k−ω and R ij −ε are compared with data from Klopman (1994).A good agreement is obtained without distinctions between the closure models.It was also of important to analyze the capacity of Code_Saturne model to reproduce the vertical profile of the Reynolds shear stress R xz = -<u'w'>.The vertical profile of R xz for the case of having only a current in the flume was tested.On Figure 4, a comparison is made between the shear stress profile obtained by Klopman (1994) and the results obtained with the k−ε, k−ω and R ij −ε models.All the turbulence closure models fit quite well Klopman's data for the whole water depth.The figures below show simulated profiles, for only a current in the channel, of the dimensionless turbulent kinetic energy (Figure 5) and the dimensionless dissipation rate (Figure 6) with the three turbulence closure models.Semi-empirical formulas (Nezu and Nakagawa 1993) were also included and used to estimate the dimensionless turbulent kinetic energy (13) and dissipation rate (14).
It can be observed that the profiles of the nondimensional turbulent intensities are quite similar to each other.The comparison of the numerical simulations with the semi-empirical curves shows, in general, the same order of magnitude, particularly close to the surface.Nevertheless, the k−ω model tends to show an increase of the values of the turbulence intensities relatively to the k−ε and R ij −ε models, mainly for the distribution of dissipation rate.The dimensionless turbulent viscosity estimated from Code_Saturne's results was also compared with experimental data obtained by Nezu and Rodi (1986) (Figure 7).
The overall features of the vertical profile of this viscosity were correctly modeled with the different turbulence closure models, but the R ij −ε model is the one that best fits the experimental data.For this reason, the latter was used to evidence the effects of the free surface boundary condition expressed by ( 11) on the ability of Code_Saturne to represent other turbulence related variables.
The zero-gradient free surface boundary condition was insufficient to get the decrease of the Reynolds shear stress towards zero (Figure 8).Only the additional empirical conditions proposed by Celik and Rodi (1984) bring the Code_Saturne to a correct behavior.The same goes for the turbulent viscosity (Figure 9).

Case with only waves, without current
After the first test with only a current in the flume, waves were generated along the positive x axis of the channel without currents.On Figure 10 we compare the mean horizontal velocity profiles in a linear scale.As it can be verified, the best agreement throughout the water depth is provided by the R ij −ε turbulence closure model .Klopman (1994) measured near the bed a change of sign on the mean horizontal velocity that becomes negative.Below this layer, where the velocity is positive, there is the wave induced streaming effect.The R ij −ε model was the only capable to reproduce this mean streaming current near the bottom, first described by Longuet-Higgins (1953) for sinusoidal water waves.As it can be seen on Figure 11, all the models were capable to reproduce quite well the distribution of the horizontal velocity amplitude profile obtained in the laboratory, not only near the bottom but also near the free surface.The k−ω model shows a little overestimation.

Waves-current interactions
When waves are superimposed on a current, significant changes are observed on the mean horizontal velocity throughout the water column.As it can be verified on Figure 12, and when comparing with Figure 3, on the upper half of the water depth the velocity shear decreases towards the free surface even reaching negative values, for the case of waves following the current.
Figure 12 shows that the R ij −ε turbulence model was the only one able to simulate the changes that occur when waves are co-flowing with the current.The first order turbulence models, k−ε and k−ω, were not able to reproduce the reduction of mean horizontal velocity values near the free surface.One possible explanation for the observed difference is the limitations of the Boussinesq assumption and a simplified vertical turbulent viscosity expression.Contrary to what happens with waves following the current, when waves counter-flow the current an increase of the velocity near the surface is observed.This behavior was observed in Klopman (1994) and others experiments as Kemp and Simons (1983).In analogous way to wavefollowing current, the R ij −ε model performed better, even if the increase of the mean velocity in the upper part of the water column was underestimated and the mean current profile in the middle of the water column was a bit overestimated (Figure 13).Nevertheless, the mean horizontal velocity profile had a good development over the water column.Moreover, the analysis of the vertical profile of the amplitude of the horizontal (orbital) velocity was made for the case of waves co-flowing (Figure 14) and counter-flowing (Figure 15) the current.A better performance for the R ij −ε approach is achieved and for the first order models there is a slight underestimation when comparing with the measurements.When waves are opposing the current (Figure 15) the improvement when using the turbulence closure model R ij −ε in Code_Saturne is clearer.

CONCLUDING REMARKS
Code_Saturne (Archambeau et al. 2004), a CFD solver based on the RANS equations, was applied to model free surface turbulent flows in a waves and current co-existing environment.There was no separation of the wave part and current part, and the two contributions are solved simultaneously.
Numerical results were compared with data from Klopman (1994), with particular attention to the vertical profiles of the mean horizontal velocities and amplitudes of the horizontal orbital velocities for three different test cases: only waves, waves following the current and waves opposing the current.A test with only a mean constant discharge was also performed, from which it was measured mean horizontal velocity and shear stress' vertical profiles A sensitivity analysis regarding the turbulence closure model in Code_Saturne appropriate for this kind of flow conditions was also made.The results were obtained without any modification of the default values of the various parameters of the selected turbulence schemes, being only imposed a boundary condition for the turbulent dissipation at the free surface.
For the case of having only a current in the channel, a more detailed study was made regarding the distribution over the water depth of some dimensionless turbulent quantities, namely, the turbulent kinetic energy, the dissipation rate and turbulent viscosity.In general, the three turbulence closure models gave quite similar results for this condition ("only current").
When superimposing the waves on the current, the second order Reynolds Stress Transport turbulence model R ij -ε SSG (Speziale et al. 1991) model has exhibited the best results, when compared with k -ω and k -ε models.The change of the vertical velocity gradient in the mean horizontal velocity profile caused by the presence of following and opposing waves in the mean flow has been well reproduced.
It was verified that the boundary condition, proposed by Celik and Rodi (1984), showed to be essential to reproduce the vertical profile of Reynolds stresses and eddy viscosity correctly.
As a general conclusion, the analysis made showed that Code_Saturne is capable of resolving the vertical structure of both mean discharge and combined flows.The change of the vertical gradient in the profile of the mean horizontal velocity caused by the presence of following or opposing waves on the mean flow has been well reproduced by the numerical model.When waves are superimposed in the same direction of the current there is a significant reduction of the mean horizontal velocity in the upper half of the water column, while when waves have opposite direction of the current the vertical shear of horizontal velocity increases.Yang et al. (2006) has stated that the cause for this behavior is due to not only the wave induced Reynolds stresses, but also to the non uniformity of the flow and secondary currents.

Figure 2 .
Figure 2. Representation of the artificial viscous beach imposed at the end of the numerical flume.

Figure 3 .
Figure 3. Vertical profiles of mean horizontal velocity for only currents.Data from Klopman (1994).

Figure 4 .
Figure 4. Vertical profiles of the Reynolds shear stress Rxz = -<u'w'> for the "only current" case using the three turbulence closure models in Code_Saturne.Data from Klopman (1994).

Figure 5 .
Figure 5. Vertical profiles of the dimensionless turbulent kinetic energy for the "only current" case using the three turbulence closure models in Code_Saturne and a semi empirical formula (13).

Figure 6 .
Figure 6.Vertical profiles of the dimensionless turbulent dissipation for the "only current" case using the three turbulence closure models in Code_Saturne and a semi empirical formula (14).

Figure 7 .
Figure 7. Vertical profiles of the dimensionless turbulent viscosity for the "only current" case using the three turbulence closure models in Code_Saturne.

Figure 8 .
Figure 8. Vertical profiles of the Reynolds shear stress Rxz = -<u'w'> for the "only current" case with and without the additional boundary condition implemented (Eq.11) and experimental data obtained by Klopman (1994).

Figure 9 .
Figure 9. Vertical profiles of the dimensionless turbulent viscosity for the "only current" case with and without the additional boundary condition implemented (Eq.11) and experimental data obtained by Nezu & Rodi (1986).

Figure 10 .
Figure 10.Vertical profiles of mean horizontal velocity for only waves.Data from Klopman (1994).

Figure 12 .
Figure 12.Vertical profiles of mean horizontal velocity for waves following current.Data from Klopman (1994).

Figure 13 .
Figure 13.Vertical profiles of mean horizontal velocity for waves opposing current.Data from Klopman (1994).

Figure 14 .
Figure 14.Vertical profiles of the amplitude of the horizontal (orbital) velocity for waves following current.Data from Klopman (1994).

Figure 15 .
Figure 15.Vertical profile of the amplitude of the horizontal (orbital) velocity for waves opposing current.Data from Klopman (1994).