# Large eddy simulation of regular waves breaking over a sloping beach

Pierre Lubin, Hubert Branger, Olivier Kimmoun

► To cite this version:

Pierre Lubin, Hubert Branger, Olivier Kimmoun. Large eddy simulation of regular waves breaking over a sloping beach. 30th International Conference on Coastal Engineering (ICCE 2006), Sep 2006, San Diego, CA, United States. pp.238-250, 10.1142/9789812709554_0021 . hal-00092703

HAL Id: hal-00092703

https://hal.science/hal-00092703v1

Submitted on 17 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# LARGE EDDY SIMULATION OF REGULAR WAVES BREAKING
OVER A SLOPING BEACH

Pierre Lubin<fnref n="1" />, Hubert Branger<fnref n="2" />, Olivier Kimmoun<fnref n="3" />

The scope of this paper is to show the results obtained for the Large Eddy Simulations (LES) of three-dimensional regular waves shoaling and breaking over a sloping beach, by solving the Navier-Stokes equations in air and water. The interface tracking is achieved by a Lax-Wendroff TVD scheme, which is able to handle with interface reconnections. The breaking processes including shoaling, overturning, splash-up and air entrainment, will be presented and discussed.

## INTRODUCTION

Great improvements have been brought to the knowledge of the hydrodynamics and the general processes occurring in the surf zone, widely affected by the breaking of the waves (Peregrine 1983, Christensen et al. 2002). Nevertheless, the turbulent flow structure is still very complicated to investigate. The interest of the numerical approach is to provide a complete and accurate description of free surface and velocity evolutions in both air and water media during the breaking of the waves, which must lead to the understanding of energy dissipation and turbulent flow structures generation processes. On the basis of the published numerical works (Lubin et al. 2006) and the experiments conducted by Kimmoun et al. (2004, 2006), a numerical study of regular waves breaking over a sloping beach is presented. Large Eddy Simulation have already been found to be a very reliable method for the unsteady simulation of the Navier-Stokes equations, as already shown in the recent works of Watanabe et al. (2005) and Christensen (2006).

The numerical methods are first introduced. We implemented up to date subgrid scale models to describe turbulence generated in both air and water media. Then, the experimental wave tank is presented, with the visualizations of the plunging breaking waves over the sloping beach. The results of two-dimensional (2D) and three-dimensional (3D) numerical simulations of the regular waves are finally compared to the experimental data. A discussion about the general behavior of the flow in both media is held. We aim at illustrating the capacity of the numerical tool to give an accurate description of the flow, in terms of instantaneous quantities, and the improvement in the numerical modelling of breaking waves since it includes the air entrainment process neglected in most previous existing models.

<footnote n="1">¹ Université Bordeaux 1, TREFLE-ENSCPB UMR CNRS 8508, 16 avenue Pey-Berland, Pessac, France</footnote>

<footnote n="2">² IRPHE UMR CNRS 6594, Aix-Marseille Université, 49 rue Joliot Curie, Marseille, France</footnote>

<footnote n="3">³ LRD Ecole Centrale de Marseille, 38 rue Joliot Curie, Marseille, France</footnote>
## MODEL AND NUMERICAL METHODS

The numerical methods have already been fully described by Lubin et al. (2006). We solve the Navier-Stokes equations in air and water, coupled with a subgrid scale turbulence model (Large Eddy Simulation – LES). The numerical tool is well suited to deal with strong interface deformations occurring during wave breaking, for example, and with turbulence modelling in the presence of a free surface in a more general way. The numerical tool has already been shown to give accurate results for coastal applications (Lubin et al. 2004). For example, comparisons with experimental and numerical results have been shown by Lubin (2004) and Helluy et al. (2005) for the practical application of the propagation and the breaking of a solitary wave over a submerged reef.

### Governing equations

On a fixed orthogonal curvilinear grid, an incompressible multiphase phase flow between non-miscible fluids can be described by the Navier-Stokes equations in their multiphase form. For each medium, a phase function $C$, or “color” function, is used to locate the different fluids standing $C = 0$ in the outer media, $C = 1$ in the considered medium. The interface between a medium and the entire domain is repaired by the discontinuity of $C$ between 0 and 1. In practice, $C = 0.5$ is used to characterize this surface. The governing equations for the Large Eddy Simulation of an incompressible fluid flow are classically derived by applying a convolution filter to the unsteady Navier-Stokes equations. The resulting set of equations reads:

$$ \nabla \cdot \mathbf{u} = 0 \qquad (1) $$

$$ \rho \left( \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} \right) = \rho g - \nabla p + \nabla \cdot \left( (\mu + \mu_T)(\nabla \mathbf{u} + \nabla' \mathbf{u}) \right) + F - \frac{\mu}{K} \mathbf{u} \qquad (2) $$

$$ \frac{\partial C}{\partial t} + \mathbf{u} \cdot \nabla C = 0 \qquad (3) $$

$$ \begin{aligned} & \rho = \rho_1 \text{ and } \mu = \mu_1 \text{ if } C \ge 0.5 \\ & \rho = \rho_2 \text{ and } \mu = \mu_2 \text{ if } C < 0.5 \end{aligned} \qquad (4) $$

with velocity $\mathbf{u}$ and pressure $p$, assuming $g$ as the gravity vector, $\rho$ as the density, $\mu$ as the viscosity, $\mu_T$ as the turbulent viscosity, $t$ as the time and $F$ as the superficial tension volume force. To deal with solid obstacles within the numerical domain, it is possible to use multi-grid domains, but it is often much simpler to implement the Brinkman theory: the numerical domain is then
considered as a unique porous medium. The permeability coefficient *K* defines the capability of a porous medium to let pass the fluids more or less freely through it. If this permeability coefficient is great (*K* → +∞), the medium is equivalent to a fluid. If it is nil, we can model an impermeable solid. A real porous medium is modelled with intermediate values of *K*. It is then possible to model moving rigid boundaries or complex geometries. To take this coefficient *K* into account in our system of equations, we thus add an extra term, called Darcy term, (μ / *K*) *u*. Model (Eqs. 1 to 4) describes the entire hydrodynamics and geometrical processes involved in the motion of multiphase media.

Large scale turbulence is described by solving the flow equations (Eqs. 1 and 2), the small scale turbulence, which is not resolved by the flow model, is taken into account through a subgrid scale model. To represent the dissipative effect of the small turbulent structures a turbulent viscosity μT is calculated with the Smagorinsky model or the Mixed Scale model. Both models are based on the resolved scales through the strain rate tensor, but the Mixed Scale model takes also into account some subgrid information. This model has been checked to be much more efficient. We also implemented a selective function to check that the velocity field is turbulent and requires a surgrid scale model to be turned on. All the details can be found in Sagaut (1998).

**Numerical methods**

A Finite-Volume method on a staggered mesh is carried out to discretize the Navier-Stokes equations and an adaptive augmented Lagrangian technique is investigated to solve the coupling between pressure and velocity in the equations of motion (Eqs. 1 and 2). The interface tracking is achieved by a Lax-Wendroff TVD scheme, which is able to handle with interface reconnections.

<figure><img src="image_2.png" /><figcaption>Figure 1. Unknown grid discretization. Black dots: pressure; black arrowheads : velocity components ; light gray dots : color function.</figcaption></figure>
A dual grid, or underlying grid (Liovic et al. 2006), is used to discretize the advection equation (Eq. 3). This technique allows to have an improved accuracy for the interface description. It can be seen in Fig. 1 that the color function is defined on a twice finer mesh than the staggered mesh used for the Navier-Stokes equations. All the references and details concerning the numerical methods can be found in Vincent et al. (2003) and Lubin (2004).

## EXPERIMENTAL CONFIGURATION

The experiments were performed in the Ecole Centrale wave tank in Marseille. A side view of the experimental apparatus is given in Fig. 2. The glass-windowed tank is 17 m long and 0.65 m wide. The water depth was set at $d = 0.735$ m. The 1/15 sloping beach was about 13 m long, starting at about 4 m away from the wavemaker. The length of the surf zone was about 3 m. Camera PIV measurements were done in twelve different locations from the incipient breaking location up to the swash zone. Fifth order Stokes waves were generated, corresponding to the analytical solution developed by Fenton (1985). The wave period was $T = 1.3$ s and wave amplitude before the sloping beach was $a = 0.07$ m. The wavelength was $L = 2.5$ m and the measured height at breaking was $H_b = 0.137$ m. The waves are observed to start breaking at about 2.65 m away from the shoreline, or 12.375 m away from the wavemaker.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Experimental wave tank.</figcaption></figure>

A sketch of a wave breaking event is displayed in Fig. 3. The wave starts breaking between windows 1 and 2 showing a short spilling phase (W2 and W3). Then a jet of liquid is rapidly ejected from the wave crest (W4) and the overturning wave front curls forward. A first splash-up is generated when the jet of liquid hits the front face of the wave (W5). We can then see a large amount of air entrained with foam and bubbles. Some other splash-ups are then generated
(W6). A roller propagates towards the shoreline, with a great air-water mixing area (W7 to W11). It can be seen in windows 7 and 8 that the bubbles are generated in the upper part of the water column, and are advected towards the bottom with a slight slanting axis (W9 and W10). The volume of entrained bubbles decreases gradually till the wave crosses the shoreline (W12), runs up before coming back. This is in agreement with the general description of Peregrine (1983), for example. More details of the experiments can be found in Kimmoun et al. (2004) and Kimmoun and Branger (2006).

<figure><img src="image_4.png" /><figcaption>Figure 3. Wave breaking in the twelve PIV interrogation windows.</figcaption></figure>

# LARGE EDDY SIMULATION OF PLUNGING BREAKING WAVES

## Wave generation procedure

Prior to the simulation of the laboratory tests, an effort has been made to implement and validate the procedure of regular and irregular wave generation developed by Lin and Liu (1999). The method consists in introducing an internal mass source function in the continuity equation for a chosen group of cells defining the source region. We adapted the analytical developments of Fenton (1989) for the fifth order Stokes wave theory, corresponding to the experimental conditions, to the source function method. In Fig. 4, we show the
velocity field around the source function, acting like a pump. Two wave trains of surface gravity waves are thus generated, as the free surface responds to a pressure increment defined within the source region cells. The two wave trains propagate in opposite directions towards the both ends of the numerical domain.

<figure><img src="image_5.png" /><figcaption>Figure 4. Velocity fields generated by the presence of the source function.</figcaption></figure>

## Computational configuration

The computational domain is 20 m long and 1.2 m high (Fig. 5). The sloping beach starts at $x = 7$ m, the source function being located at $x = 3.5$ m and $z = 0.3675$ m (the center of the source region is at $d/2$). The numerical beach is considered as an impermeable solid obstacle, the permeability coefficient $K$ being initialized at zero (Eq. 4). The source region is 0.06 m wide and 0.0735 m high. The area and the location of the source function have been designed applying the rules described by Lin and Liu (1999). 520 000 mesh grid points are used to discretize the numerical domain (2000 x 260). The grid system is uniform along the wave propagation direction, with a constant cell length $\Delta x = 1.10^{-2}$ m, and non-uniform in the vertical direction, the minimum cell height being $\Delta z = 2.5 \cdot 10^{-3}$ m. These values have to be divided by two, in both directions, for the free surface description thanks to the dual grid. The time step is chosen to ensure a Courant-Friedrichs-Levy condition less than 1, necessary for the explicit advection of the free surface.

As already said, two wave trains are generated and propagate in opposite directions towards the both ends of the numerical domain. An open boundary condition is thus set at the left side of the numerical domain to let the outgoing wave exit the numerical domain. In order to ensure that no numerical reflection occurs at the left side of the numerical domain, a sponge layer is set in addition to the open boundary condition. It consists in a region where the permeability coefficient $K$ is chosen such that the outgoing wave train is properly attenuated before reaching the open boundary. The simulations have been carried out on
fast PC's. About twenty breaking waves have been simulated. Only
instantaneous quantities are presented and discussed.

<figure><img src="image_6.png" /><figcaption>Figure 5. Sketch of the numerical configuration for the regular wave breaking simulation. Dashed line: mean water depth d = 0.735 m. The toe of the 1/15 sloping beach starts at x = 7 m.</figcaption></figure>

## 2D numerical results

In Fig. 6, we present the numerical results corresponding to the 10th breaking wave. The figures correspond to the twelve experimental pictures taken at the different locations in the surf zone (Figs. 3). A general good agreement has been found in terms of wave propagation. The main observed difference is that we miss the first short spilling event.

Some elements can be brought to discuss this discrepancy. On the basis of the experimental measured waveheight at breaking $H_b$, the surf similarity parameter has been evaluated to be $\xi = 0.285$, which usually characterizes a spilling breaker (Battjes 1988). Spilling breaking waves involve a complex combination between the wave propagation and the vorticity generation based on perturbations appearing on the crest of the steepening face of the wave. Moreover, the measured height of the white cap is about 1 mm (Fig. 3 – windows 2 and 3). So, numerically, the mesh grid refinement we use is obviously not sufficient enough to be able to capture this feature. The numerical diffusion is also probably responsible, even if it has been checked to be very low (the numerical diffusion of the free-surface is contained in three mesh grid cells).

Once the front face of the crest steepens and becomes vertical (Fig. 6 - windows 1 to 4), a jet of water is projected (Fig. 6 - between windows 4 and 5), and free falls down forward into a characteristic overturning motion. The plunging breaking wave is then responsible for the generation of jet-splash
cycles, this, in turn, being responsible for the generation of a sequence of large scale coherent vortices.

<figure><img src="image_7.png" /><figcaption>Figure 6. Numerical results for a wave breaking at the twelve PIV interrogation windows locations (10th wave).</figcaption></figure>

A first high splash-up is generated (Fig. 3, window 6 - Fig. 6 window 7) and the successive splash-ups are also simulated. We have air entrapped and entrained in the water. And once the wave moves into a roller propagating towards the shoreline (Fig. 6 - windows 8 to 12), we have very similar results for the free surface description. It is distorted and very dynamic (Fig. 3 - windows 8 to 12). Nevertheless, some discrepancies appear. Due to the coarse mesh grid resolution, the dislocation of the gas pockets into small bubbles cannot be simulated, even if, in the numerical results, the gas pockets correspond to some air-water mixing zones observed in the experimental pictures. Turbulence is associated with air entrainment, which is responsible for wave energy damping in the surf zone. In the experiments, it appears that the entrained air bubbles are contained mostly in the large structures and diffused towards the bottom due to the eddies. The rate of energy dissipation is increased with the bubble penetration depth and strong vertical motions are induced by the rising air bubbles. These mechanisms are mostly 3D, which cannot be taken into account in a 2D numerical simulation.

In Fig. 7, we present the numerical velocity fields corresponding to the 10th breaking wave. The figures correspond to three experimental pictures taken at the different locations in the surf zone.
<figure><img src="image_8.png" /><figcaption>Figure 7. Comparison between instantaneous PIV and numerical velocity fields at different locations.</figcaption></figure>

A high velocity region is located at the crest of the wave. A high splash-up is rising, with high velocities directed upward and towards the wave propagation direction, which generates counter-rotative vortices, as observed by Bonmarin (1989). This process is observed to repeat, each successive splash-up being weaker than the preceding one. Some large volumes of air and water are put into rotation. Co-rotative vortices are then generated by these successive splash-ups, as observed by Miller (1976) or Sakaï et al. (1986). High velocities are located near the free-surface, due to the jet-splash cycles, during the bore propagation.

Considering these simulations as a validation step, our numerical model gives very satisfactory and encouraging results for this 2D configuration. Nevertheless, wave breaking is a 3D two-phase turbulent problem, so the first numerical results will be presented hereafter, consisting in a first attempt.
3D numerical results

<figure><img src="image_9.png" /><figcaption>Figure 8. Numerical results for 3D waves breaking (14th wave). Color function $C = 0.5$.</figcaption></figure>

The same initial condition than previously detailed is used, spread along the transverse direction. The numerical wave tank is 30 cm wide. But the numerical domain has been shortened and is discretized with 660 000 mesh grid points (550 x 60 x 20) to be still able to run the simulation on a fast PC with a reasonable accuracy. The numerical domain is now 15.5 m long. The grid system is uniform along the wave propagation direction, with a constant cell length $\Delta x = 2.8 \cdot 10^{-2}$ m, and non-uniform in the vertical direction, the minimum cell height being $\Delta z = 8.4 \cdot 10^{-3}$ m. A constant cell width $\Delta y = 1.5 \cdot 10^{-2}$ m is used.
These values have to be divided by two, for all directions, for the free surface description thanks to the dual grid.

<figure><img src="image_10.png" /><figcaption>Figure 9. Numerical results for 3D waves breaking (end of the 14th wave). Color function C = 0.5.</figcaption></figure>

About fifteen breaking waves have been simulated. The number of grid points is obviously not sufficient enough again to take into account the dislocation of the pockets of gas into small bubbles, but the flow dynamics is
correctly described in the surf zone with air entrainment. A plunging breaker type is still observed. Figs. 8 and 9 present the numerical results corresponding to the 14th wave breaking. The amount of air entrained in the water is less important than previously observed in the 2D numerical results due to the less accurate description of the plunging jet, responsible for a great quantity of air to be entrapped.

## CONCLUSION

The numerical results presented in this paper concerns instantaneous quantities, simulating 2D, and then 3D, regular waves breaking over a sloping beach. Numerical simulations are now undertaken to calculate phase averaged and fluctuating quantities. Turbulence will then be studied quantitatively and compared with the experiments. Our model was found to be reliable to describe correctly and accurately the complicated two-phase flow interactions that happens when waves break. The breaking process, in terms of wave overturning and splash-up occurrence, is in accordance with the general observations given in the literature. The air entrainment is described, which is important as it plays a great role in the energy dissipation process.

But 3D numerical simulations remain very demanding and time consuming. Two-phase flows simulations and turbulence modelling require fine mesh cells to be very accurate, even if Large Eddy Simulation is supposed to allow us to save some mesh points. Some theoretical work is also led to model what happens in the vicinity of the color function $C = 0.5$ (Lubin 2004, Labourasse et al. 2006).

## ACKNOWLEDGMENTS

The authors would like to acknowledge the financial and scientific support of the French INSU - CNRS (Institut National des Sciences de l'Univers - Centre National de la Recherche Scientifique) program IDAO (« Interactions et Dynamique de l'Atmosphère et de l'Océan » and the IDRIS - CNRS (« Institut du Développement et des Ressources en Informatique Scientifique ») and CINES (« Centre Informatique National de l'Enseignement Supérieur ») for the material support (project ter 2237).

## REFERENCES

- Battjes, A. 1988. Surf-zone dynamics. *Ann. Rev. Fluid Mech.*, 20, 257-293.

- Bonmarin, P. 1989. Geometric properties of deep-water breaking waves, *J. Fluid Mech.*, 209, 405-433.

- Christensen, E. D., D. J. Walstra, and N. Emarat. 2002. Vertical variation of the flow across the surf zone, *Coastal Engineering*, 45, 169-198.

- Christensen, E. D. 2006. Large eddy simulation of spilling and plunging breakers, *Coastal Engineering*, 53, 463-485.

- Fenton, J. D. 1985. A fifth-order Stokes theory for steady waves. *Journal of Waterway, Port, Coastal and Ocean Engineering*, 111, 216-234.
Helluy, P., F. Gollay, S. T. Grilli, N. Seguin, P. Lubin, J.-P. Caltagirone, S. Vincent, D. Drevard and R. Mercer. 2005. Numerical simulations of wave breaking. *Mathematical Modelling and Numerical Analysis Mathematical Modelling and Numerical Analysis*, 39 (3), 591-608.

Kimmoun O., H. Branger and B. Zucchini. 2004. Laboratory PIV measurements of wave breaking on a beach. *Proceedings of 14th International Offshore and Polar Engineering Conference*, 3, 293-298.

Kimmoun K. and H. Branger. 2006. A PIV investigation on laboratory surf-zone breaking waves over a sloping beach, *corrected and submitted to J. Fluid Mech*.

Labourasse, E., D. Lacanette, A. Toutant, P. Lubin, S. Vincent, O. Lebaigue, J.-P. Caltagirone, P. Sagaut. 2006. Towards Large Eddy Simulation of isothermal two-phase flows: governing equations and a priori tests, *under press for International Journal of Multiphase Flow*.

Lin P. and P. L.-F. Liu. 1999. Internal wave-maker for Navier-Stokes equations models, *Journal of Waterway, Port, Coastal and Ocean Engineering*, 125, 322-330.

Liovic P., M. Rudmanb, J.-L. Liowc, D. Lakehal and D. Kothe. 2006. A 3D unsplit-advection volume tracking algorithm with planarity-preserving interface reconstruction, *Computers and Fluids*, 35 (10), 1011-1032.

Lubin, P. 2004. *Large eddy simulation of plunging breaking waves*, PhD thesis (in English), University of Bordeaux.

Lubin, P., S. Vincent, J.-P. Caltagirone and S. Abadie. 2004. Three-dimensional large eddy simulation of vortices induced by plunging breaking waves, *Proceedings of 29th International Conference on Coastal Engineering*, ASCE, 344-356.

Lubin, P., S. Vincent, S. Abadie and J.-P. Caltagirone. 2006. Three-dimensional Large Eddy Simulation of air entrainment under plunging breaking waves, *Coastal Engineering*, 53, 631-655.

Miller, R. L. 1976. Role of vortices in surf zone prediction: sedimentation and wave forces, Soc. Econ. Paleontol. Mineral. Spec. Publ., 24, Ed. R. A. Davis, R. L. Ethington, 92-114.

Peregrine, D. H. 1983. Breaking waves on beaches, *Ann. Rev. Fluid Mech.*, 15, 149-178.

Sagaut, P. 1998. Large eddy simulation for Incompressible Flows, Springer, Verlag.

Sakaï, T., T. Mizutani, H. Tanaka and Y. Tada. 1986. Vortex formation in plunging breaker, *Proceedings of 20th International Conference on Coastal Engineering*, ASCE, 711-723.

Vincent, S., J.-P. Caltagirone, P. Lubin and T. N. Randrianarivelo. 2003. An adaptive augmented Lagrangian method for three-dimensional multi-material flows, *Computers and Fluids*, 33, 1273-1289.

Watanabe, Y., H. Saeki and R. J. Hosking. 2005. Three-dimensional vortex structures under breaking waves, *J. Fluid Mech.*, 545, 291-328.