THREE-DIMENSIONAL NUMERICAL SIMULATIONS OF AERATED VORTEX FILAMENTS UNDER PLUNGING BREAKING WAVES

The scope of this work is to present and discuss the results obtained from simulating three-dimensional plunging breaking waves. A numerical tool based on the Navier-Stokes equations is used to describe the plunging breaking processes including overturning, splash-up and the occurrence of air entrainment. Initial 3D conditions corresponding to unstable periodic sinusoidal waves of large amplitudes in periodic domains are used to study the air entrainment occurring when waves break. The numerical results highlight the major role of this phenomenon in the energy dissipation process, through a high level of turbulence generation. The numerical model represents a substantial improvement in the numerical modelling of breaking waves since it includes the air entrainment process neglected in most previous existing models.


INTRODUCTION
In the last three decades, a large amount of studies have been devoted to improving the description of the hydrodynamics and the general processes occurring in the surf zone, widely affected by the breaking of the waves.
Depending on the type of breaker, three types of large-scale coherent structures can be found in breaking waves (Zhang & Sunamura 1994).Occurring several times in a single plunging breaker, the jet-splash cycle, are responsible for the generation of a sequence of large-scale horizontal vortices.These structures, having a horizontal axis of rotation, have been described by many observations (Sawaragi & Iwata 1994;Miller 1976;Bonmarin 1989;Kimmoun & Branger 2007).For the spilling breaker condition, Nadaoka et al. (1989) investigated the flow field under the turbulent bore propagating towards the shoreline.Behind the wave crest, the flow structure changes rapidly into obliquely downward stretched three-dimensional eddies (obliquely descending eddies).Kubo & Sunamura (2001) revealed recently that a new type of large-scale turbulence, named the "downburst", can be found in the breaker zone aside from the oblique vortices.Ting (2006Ting ( , 2008) ) also identified these downbursts of turbulence descending from the free-surface.
Few numerical works have been dedicated to reproducing these three-dimensional large features occurring when waves break (Watanabe et al. 2005;Lubin et al. 2006;Iafrati 2009;Lakehal & Liovic 2011).Moreover, only a limited number of researchers have also focused attention on some smallerscale processes (Longuet-Higgins 1995).Narayanaswamy & Dalrymple (2002) experimentally presented evidence of "fingers" appearing at the tip of the plunging jet, prior to impact.More recently, Saruwatari et al. (2009) studied numerically the formation of fingers and scars at the surface of secondary planar jets and suggested that, as the influence of surface tension increases, the jet surface is prevented from being scarified and fingered.Handler et al. (2012) experimentally investigated the generation of coherent elongated structures behind breaking waves.
On 6 th May 2009, BBC Worldwide released a short clip of a large breaking wave filmed in slow motion (BBC 2009).The beautiful breaking wave, filmed from underwater, revealed for the first time three-dimensional coherent structures: tornado-like vortical tubes, connecting the splash-up and the main tube of air.
However, the wave breaking phenomenon remains a very challenging fluid mechanics problem, due to turbulence and aeration interactions.Performing numerical simulations of three-dimensional breaking waves requires a large number of mesh grid nodes, robust and accurate numerical methods, and long CPU time calculations to compute the hydrodynamics from the largest to the smallest length and time scales (Lubin et al. 2011).Recent progress in computational capacities allowed us to run fine three-dimensional simulations giving us the opportunity to observe for the first time these fine vortex filaments generated during the early stage of the plunging wave breaking process, and the subsequent air entrainment (Lubin & Glockner 2013).The scope of this paper is to present and discuss the visualizations of these coherent eddy structures.The numerical results will be detailed to educe and explain the formation and evolutionary dynamics of the vortex filaments, and to explore the role of the eddies in air entrainment.

MODEL AND NUMERICAL METHODS
We solve the Navier-Stokes equations in air and water, coupled with a subgrid scale turbulence model.The numerical tool is well suited to deal with strong interface deformations occurring during wave breaking, for example, and with turbulence modeling in the presence of a free-surface in a more general way.

Governing equations
An incompressible multiphase phase flow between non-miscible fluids can be described by the Navier-Stokes equations in their multiphase form.The governing equations for the Large Eddy Simulation (LES) of an incompressible fluid flow are classically derived by applying a convolution filter to the unsteady Navier-Stokes equations.In the single fluid formulation of the problem, a phase function C, or "color" function, is used to locate the different fluids standing C = 0 in the outer medium, C = 1 in the considered medium.The interface between two media is repaired by the discontinuity of C between 0 and 1.In practice, C = 0.5 is used to characterize this surface.
The resulting set of equations reads (Eqs.1-4): (1) with velocity u and pressure p, assuming g as the gravity vector, ρ as the density, µ as the viscosity, µ T as the turbulent viscosity and t as the time.
The magnitudes of the physical characteristics of the fluids are defined according to C in a continuous manner as: (1 ) where ρ 0 , ρ 1 , µ 0 and µ 1 are the densities and viscosities of fluid 0 and 1 respectively.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 Mixed Scale model (Sagaut 1998).Model (Eqs. 1 to 4) describes the entire hydrodynamics and geometrical processes involved in the motion of multiphase media.

Numerical methods
Time discretization of the momentum equation is implicit and an Euler scheme is used.The velocity/pressure coupling under the incompressible flow constraint is solved thanks to the time splitting pressure correction method (Goda, 1979).The equations are discretized on a staggered grid by means of the finite volume method.The space derivatives of the inertial term are discretized by a hybrid upwind-centered scheme.The viscous term is approximated by a second order centered scheme.
The interface tracking is achieved by a Volume Of Fluid method (VOF) and a Piecewise Linear Interface Calculation (PLIC) (Youngs, 1982;Scardovelli & Zaleski, 1999).This method has the advantage to build a sharp interface between air and water.Since the color function is not defined on each point where viscosities and densities are needed for the Navier-Stokes discretization, the physical characteristics are interpolated on the staggered grid.We use a linear interpolation to calculate the density on the velocity nodes, whereas a harmonic interpolation is used for the viscosity.
The MPI library is used to parallelize the code.The mesh is partitioned into equal size subdomains to ensure load balancing.Communications between processors are also minimized (Ahusborde & Glockner, 2011).The HYPRE parallel solver and preconditioner library is used to solve the linear systems (Falgout et al., 2006).The prediction and correction steps are solved, respectively, thanks to a BiCGStab solver, associated with a point Jacobi preconditioner, and a GMRES solver, associated with a multigrid preconditioner.The numerical code has already been extensively verified and validated through numerous test-cases including mesh refinement analysis (Lubin et al., 2006(Lubin et al., , 2010(Lubin et al., , 2011;;Poux et al., 2011).The accuracy of the numerical schemes and the conservation laws of mass and energy in the computational domain have been accurately verified.

Initial and boundary conditions
We use initial conditions corresponding to a single unstable periodic sinusoidal wave of large amplitude.This somewhat artificial wave breaking configuration has already been documented in several previous studies and proved to be effective at simulating all types of breaking waves (Abadie et al. 1998;Chen et al. 1999;Lubin et al. 2006;Iafrati 2009;Hu et al. 2012).
The rectangular calculation domain is periodic in the wave propagation direction and one wavelength long (Figure 1).A free slip boundary condition is imposed at the lower limit, and a free boundary condition in the upper limit.The reference variables of the initial incident wave are the celerity c (m.s -1 ), the period T (s), the wavelength L (m), the water depth d (m), the waveheight H (m) and the densities and viscosities of air and water (kg.m -3 and kg.m -1 .s - ).The flow motion is characterized by the Reynolds number, Re = ρ w cL/μ w , the density ratio, ρ a /ρ w , the viscosity ratio, μa/μw, the initial steepness, H/L, and the dispersion parameter, d/L.In our specific case, the overturning motion is controlled only by these last two initial parameters, H/L and d/L, which are chosen such as the initial wave cannot remain steady as the initial velocity field in water is not in equilibrium with the initial wave profile.H/L mainly controls the instability growth speed, so the time separating the start of the computation and the breaking.It makes this approach interesting as we are then able to study any breaker type by varying only these flow parameters (Lubin et al. 2006;Iafrati 2009).It is also a convenient configuration as we can study accurately the breaking phenomenon in a smaller numerical domain, compared with simulations of shoaling waves breaking over a sloping beach which are very demanding in terms of grid mesh points (Lubin et al. 2011).The purpose of the paper work is to simulate 3D plunging breaking waves and investigate the generation of the vortex filaments in the resulting flow.The three-dimensional numerical domain is discretised into 1024x400x200 non-regular Cartesian cells (almost 82 millions of mesh grid cells), and partitioned into 576 subdomains (one processor per subdomain).The grid is evenly distributed in the longitudinal and transverse directions, giving a mesh grid resolution of Δx = Δy = 1.10 -4 m.In the vertical direction, the grid is clustered near the interface to resolve the small structures of the flow, with a constant grid size Δz = 1.10 -4 m below d+H=2.A non-regular grid resolution is used in the air with a maximum mesh grid size at the top of the numerical domain.The real air and water physical properties are used.At the initial time of the simulation, the water velocity field in the wave is obtained from the linear theory and the air is at rest.Hydrostatic pressure is initialized in the whole numerical domain.The computing time was approximately 48 hours for a simulated time of 0.88 ms.The time step was chosen to ensure a Courant-Friedrichs-Levy number less than 0.3.

Results and discussion
On the basis of the previous works (Lubin et al. 2006;Lubin & Glockner 2013), we simulated plunging breaking waves, ranging from weak to energetic events.According to previous numerical simulations (Lubin et al. 2006), for a given dispersion parameter value, the greater the wave steepness is, the stronger the plunging breaking event will be.The values used for the numerical simulations are presented in Table 1.We choose to identify and educe the vortex filaments motions and their evolutions in space and time using the Qcriterion.The vortex filaments are found under plunging breaking waves (Table 2), while no vortex filament is observed under the weakest plunging breaking wave.From all the simulations performed, it has been found that the vortices are generated when a plunging breaking event is observed with a plunging jet impacting far from the incident wave crest.

Yes
As many authors indicated, the plunging jet does not penetrate the forward face of the wave, whatever the position of the plunge point or the angle between the falling crest and the front face of the wave (Peregrine 1983;Jansen 1986;Lin & Hwung 1992;Kiger & Duncan 2012).The plunge point is shown to be located above the still water level.The tongue of water first hits the front face of the breaking wave.From this point, the tongue of water separates, due to the low speed of penetration compared to the flow velocity in the jet.One part of the liquid forms the upper part of the splash-up, while the other goes backwards around the main pocket of entrapped air.The initial jet of water is almost totally reflected when it hits the front face of the breaking wave (Abadie et al. 1998;Yasuda et al. 1999;Lubin et al. 2006;Dalrymple & Rogers 2006).
When the jet impinges on the face of the wave, it eventually creates a line of craters upstream from the plunge point, which will further develop into the splash-up (Fig. 2 a).The craters do not penetrate deeply, but expand due to the adverse flow from the forward face.These craters were observed to develop and deform, generating aerated cavities.The craters are deformed as they grow due to the increasing pressure between the colliding jet and the face of the wave.The bottoms of the craters are then pulled downward.A swirl then develops due to the combined effect of the downward flow of the water inside the walls of the craters and the shear created by the forward motion of the jet and the backward motion of the face of the wave.A very strong spiraling tendency is then generated and, eventually, there is enough torque to create a vortex.The shapes of the craters are non-uniform, so there is no preferred direction of rotation.The water from the impinging jet closes over the crater, entrapping a pocket of air, and separates to form the inner wall of the main tube of air.When the splash projects forward from the impact point, the notch creates an angle responsible for the tilting of the nascent vortex filaments (Fig. 2 b).As the liquid in the jet clearly separates, the vortex filaments are elongated, the upstream end of the vortex filaments remain attached to their generation points at the upper surface of the developing splash-up.Once the spiral motion has started, it tends to entrain and merge with the surrounding air pockets and bubbles.Figure 2 c shows that the upstream end of the vortex, or its ``mouth``, is opened at the extremity of the developing splash-up.The downstream end is stretched, entrained towards the horizontal tube of air and wraps around it (Fig. 3).
The eddies at the early stage of the breaking are disorganized, whereas during the splash-up development phase a more regular pattern can be observed with quasi-streamwise coherent vortices (Fig. 4).The vortex filaments now develop into 3D streamwise coherent structures.The upstream ends remain attached to the generation points at the upper surface of the splash-ups, while the downstream ends are connected to the main tube of air.The rotating structures are thus clearly identified as soon as the jet enters the forward face of the wave.The vortex filaments are found to be located at the exact boundary between the water from the impacting jet and the water from the forward face of the wave, in the strain region.The complete description of the generation mechanism is provided by Lubin & Glockner (2014).
We observed on the back of the plunging jets some local surface patterns (Figure 2), which are expected to be linked to the striations or the fingers detailed respectively by Longuet-Higgins (1995) and Narayanaswamy & Dalrymple (2002).This could be a key phenomenon explaining the transverse distribution of the vortex filaments.Some existing analogous streamwise vortical structures can be observed in shear flows.The main difference is that we identified and explained how these vortex filaments were individually generated in the chaotic impact of a mass of water.The structures are observed to be aerated, as their mouths are attached to the free-surface while the other end is linked to the main tube of air.We do not noticed any bending at their extremities, nor continuous linking to neighbouring filaments to form loops (Watanabe et al. 2005;Saruwatari et al. 2009).The common point is the shear region in the saddlepoint, responsible for the stretch-and-intensification process of the filaments.Each extremities of the vortex filaments have been seen to be connected: one end is attached to the main tube of air and wraps around it, and the other end is opened at the free-surface, at the upper part of the splash-up.Moreover, aeration of the vortex cores make the filaments hollow.

CONCLUSIONS
We observed that these vortex filaments initiated when the tip of the jet touched down the forward face of the wave.Then, the liquid in the jet separated, one part of the liquid formed the upper part of the splash-up, the other went around the main entrapped tube of air.The vortex filaments were then elongated from their generation points, stretched from the developing splash-up towards the other ends of the filaments, bending and wrapping around the main aerated cavity.Some of the vortex filaments have been observed to entrain air, while some structures interacted with each other, forming larger vortex filaments and sometimes coiling.
The streamwise vortex filaments have been detected under plunging breaking waves only, with an upstream obliqueness of about 50°.They have been shown to be independent of the breaker intensity.Their lateral distribution at the generation seemed to be driven by the striations found on the back of the plunging jet.
Regular waves should be studied in future works.We simulated a single breaking event in a periodic domain, while Watanabe et al. (2005) simulated regular waves breaking one after another.So, we have to evaluate how each wave would interacts with the remaining vorticity left by the previous breaking wave.The overall vorticity field should be more complex.For this reason, we also limited our discussion to the generation of the vortex filaments and chose not to develop their interactions after wrapping around the main spanwise tube of air.Future works will be dedicated to address the generation of the other types of 3D structures encountered under breaking waves (obliquely descending eddies and downbursts).The generation of the striations on the back of the plunging jet is also a subject of future interest.

Figure 1 .
Figure 1.Sketch of the initial condition for the unstable periodic sinusoidal wave of large amplitude.The interface is located at C=0.5.Note that the Navier-Stokes equations are solved in air and water.H/L=0.13,d/L=0.13,t=0 s.

Figure 2 .
Figure 2. Sequence of pictures presenting the evolution of a single vortex filament, from the touchdown of the plunger to the beginning of the splash-up generation, at times t = 0.1499 s, 0.1587 s and 0.1698 s.The vortex filament is isolated in a 3D strip extracted from the whole numerical domain for a clearer identification of the structure.

Figure 3 .
Figure 3. Sequence of pictures presenting the evolution of the vortex filaments, from the touchdown of the plunger to the beginning of the splash-up generation, for H/L=0.13 and d/L=0.13.The vortex filaments are isolated in a 3D strip extracted from the whole numerical domain for a clearer identification of the structures.The free-surface is identified with the isocontour of the phase function C=0.5 (in blue).The isocontour of the axial vorticity is shown in yellow.

Figure 4 .
Figure 4. Coherent vortical structures underneath the plunging breaking wave (isocontour of the phase function in blue, isocontour of the Q-criterion in green).