Wave-breaking and generic singularities of nonlinear hyperbolic equations
aa r X i v : . [ phy s i c s . f l u - dyn ] A p r Wave-breaking and generic singularities of nonlinear hyperbolic equations
Yves Pomeau
Los Alamos National Lab, CNLS, Los Alamos, NM 87545 USA.
Martine Le Berre
Laboratoire de Photophysique Mol´eculaire, Bat.210, 91405 Orsay, France.
Philippe Guyenne
Department of Mathematical Sciences, University of Delaware, Newark DE 19716-2553, USA.
Stephan Grilli
Department of Ocean Engineering, University of Rhode Island, Narragansett, RI 02882, USA. (Dated: October 22, 2018)
Abstract
Wave-breaking is studied analytically first and the results are compared with accuratenumerical simulations of 3D wave-breaking. We focus on the time dependence of various quantitiesbecoming singular at the onset of breaking. The power laws derived from general arguments and thesingular behavior of solutions of nonlinear hyperbolic differential equations are in excellent agreementwith the numerical results. This shows the power of the analysis by methods using generic conceptsof nonlinear science.
PACS numbers:
I. INTRODUCTION
Among the myriad of natural phenomena dominated by nonlinearity, wave-breaking is quite special. It is ubiquitousalong most oceanic shorelines, and although progress has recently been made in its understanding, its full explanationremains somewhat shrouded in mystery. This makes wave-breaking an outstanding topic in the study of nonlinearphenomena and, we believe, one that is particularly well suited to this “Open problems section of Nonlinearity”.Wave-breaking also plays a significant role in the larger scale dynamics of the Oceans [1], particularly with respectto their interaction with the atmosphere. For instance, the swell provides nucleation seeds for the rain and windmomentum exchange with water, mostly through nonlinear processes. Among such processes, wave-breaking couldalso be significant, if not dominant, to limit the growth of the instability responsible of the formation of waves.The linear theory of waves has a long and interesting history, going back to the Principia. In book 2 Newton showsthat the velocity of water waves in the infinite depth limit is proportional to √ gλ , with g the acceleration of gravityand λ the wavelength. Years later, Bernoulli derived the shallow water wave velocity, applicable to waves with awavelength much larger (typically more than 20 times) than the water depth. The interpolation between these twolimits is due to Lagrange. Later, in a very important paper, Stokes defined the range of applicability of the lineartheory in the shallow water limit. Thereafter Riemann’s method of solution of quasilinear equations showed theoccurrence of a finite time singularity for simple waves in the shallow water case. Among the open problems posedby the occurrence of such a singularity, we shortly discuss below the two existing methods proposed to regularize thewave profile to avoid the occurrence of an overhang. These methods introduce the next order space derivative in theshallow water equations (i.e., so-called dispersive effects) and thus preclude wave-breaking from occurring.The relationship between wave-breaking and the nonlinear structure of the governing equations is reconsideredbelow, where we show that strong nonlinearities dominate the wave dynamics, both in the transverse and longitudinaldirections. Starting from generic equations, whose solution can develop discontinuities (inviscid Burgers, or pressure-less fluid equations), we first show in a one-dimensional space (1D), that a very simple (once known as [2]) solution u ( x, t ) of this partial differential equation does yield a wave-breaking (i.e., overturning) process. Defining u as the freesurface elevation, we study its dynamical behavior and, more specifically, the growth of the overhang domain withtime.Extending the pressure-less fluid equation to a two-dimensional space (2D), we obtain a generic equation for u ( x, y, t )showing how the “overhanging domain” of the overturning wave extends in both spatial directions, parallel and normalto the wave. In that case, the boundary of the overhanging region, referred to as the “apparent contour” (i.e., theset of points where the tangent to the water surface is vertical), becomes a skewed curve, when drawn on the watersurface. We find that this curve spreads as t / in the direction of wave propagation and as t / along the wavecrest ( t = 0 being the instant where wave overturning is initiated, yielding an overhanging region). These resultswill be shown to be in good agreement with those of numerical simulations of a 3D overturning wave over a slopingunderwater ridge, based on the solution of Euler equations with free surface boundary conditions [6]. Moreover thelateral spreading as t / was recently found experimentally[7] for waves generated by a wave-maker, and progressingon an horizontal bottom.These time evolution laws are established first for a pressure-less fluid and then are extended in section III to theshallow water case. The scaling laws for wave-breaking, in the directions perpendicular and parallel to the crest, willbe derived from solutions of fluid-mechanical equations. It is true that this wave-breaking process could be seen asan example of a metamorphosis of a surface studied in catastrophe theory, named a ’sickle’ or ’lips’. In particular,Arnold [4], referring to physical situations completely different of ours, notices that the width of a ’sickle’ should growgenerically like ( t − t ) / , t being the time of the first singularity. This law of growth is the same as the one wefind for the lateral spreading of the overturning domain. However catastrophe theory is very general and cannot beconsidered as providing any kind of mathematical method for solving the equations of fluid mechanics. Therefore wedid not feel useful to use this theory explicitly in the present work and we relied instead on explicit solutions of thefluid equations. Moreover, to the best of our knowledge, catastrophe theory has never been used in investigations ofthe wave-breaking problem. This being a contribution to the Open problems section of the Journal, in section 4 wecomment on various aspects of the problem that do wander off the main track. For instance, we suggest an extensionof the boundary integral method (e.g., such as used in [6]) to deal with the wave-breaking problem, beyond the timewhen the plunging jet impinges on the free surface underneath, and briefly discuss some features of what we call theHertz-fluid contact problem. II. FROM 1D TO 2D WAVE-BREAKING
In this section, we first introduce a generic singular solution to a nonlinear hyperbolic partial differential equation(PDE), in a one-dimensional space (1D). We then show that such solutions of the PDE can be represented by a curvethat is continuously deforming and, hence, can be used as a simple geometric model for predicting the occurrence ofan overhanging region, similar to the plunging jet in an overturning (i.e., breaking) water wave. We finally extendthis equation to a two-dimensional space (2D), which is a more physically relevant case.
A. Wave-breaking in one dimension
Consider the equation ∂u∂t + u ∂u∂x = 0 . (1)It has [2] the implicit solution u ( x, t ) = u ( x − ut ) , (2)where u ( x ) is the initial condition. Let f ( . ) be the inverse function of u ( . ) (all functions assumed smooth). Therefore: f ( u ( x, t )) = x − ut. (3)This solution u ( x, t ), as given by equation (3), becomes multivalued whenever ∂x∂u = 0, the derivative being takenat constant t . From equation (3), this yields t + d f d u = 0 . (4)The smallest t when this has a root, for a given f ( u ), is when this root is stationary with respect to t , that is whend t = 0, with t and u related by (4). This yields d t = − d u d f d u = 0, or simply d f d u = 0. Therefore, at the singularity, f has a vanishing second derivative, a standard result. Furthermore it can also be assumed that the Taylor seriesexpansion (TSE) of f has no linear term (amounting to take u = 0 as the value of u where f ( . ) has an inflectionpoint), since a constant u can be absorbed into a redefinition of the time origin in equation (3). Let us then define t = 0 as the first time when the solution becomes singular, to avoid having to write ( t − t ∗ ) / etc. instead of t / ,with t ∗ time of the singularity. Near the singularity (assumed to occur at t = 0) one can thus expand (after convenientscaling) f ( u ) like f ( u ) = − u + au + ... , with a finite constant. Once put into equation (3) this gives to the 4th-order: au = u + x − ut. (5)For short times , Equation (5) implies that u scales like t / and x like t / , which makes all terms on the right-handside of Equation (5) of the same order, t / . The left-hand side term is - with the same scaling law - of order t andtherefore negligible compared to the right-hand side; this applies as well to higher-order terms in the TSE of f ( u )near the inflection point. Assuming that the beginning of an overhanging region can be described by the lower-orderterms, the ’universal’ equation describing ’wave-breaking’ reads:0 = u + x − ut. (6)Note that, in this derivation, the coefficient of u of the TSE of f ( u ) near zero is negative and thus a singularity willappear for positive times. For negative times, Equation (6) only has one real root while it has three such roots forpositive times.Therefore the curve u ( x, t ) is everywhere single-valued if its slope is negative at x = 0, which happens for t negative.By contrast, u ( x, t ) is triple-valued in a range of values of x , if t is positive. This solution can be put in a self-similarform by eliminating one variable. The scaling laws for the transition from single to triple-valued results are u ∼ p | t | and x ∼ | t | / . Therefore, the ’natural’ solution of Equation (5) is of the form u ( x, t ) = p | t | U ( ζ ), with ζ = x | t | − / and U a numerical function. Nevertheless, this is not sufficient, because the function U ( . ) takes different valuesdepending on whether t is negative, zero, or positive. This can be taken into account formally by introducing adiscrete argument in the function U ( . ), namely an index i with three possible values −
1, 0 and 1, such that i = − t is negative, i = 0 if t = 0 and i = 1 if t is positive. This defines three real numerical functions U i ( ζ ), such that U − ( ζ ) is the real root of 0 = U − + ζ + U − , and U = − ζ / . The multivalued function U is any real root of0 = U + ζ − U , with three possible values in the range − √ < ζ < √ and one otherwise.This completely defines the self-similar solution. B. Wave-breaking in more than one space dimension
Our approach to the formation of singularities in solutions of hyperbolic equations can be extended to nonlineardifferential equations with more than one spatial variable. A simple extension of the nonlinear equation (1) provides amodel for the occurrence of singularities for time-dependent functions of two coordinates of space ( x, y ). Let u ( x, y, t )and v ( x, y, t ) be the two Cartesian horizontal components of the velocity field of a fluid, with ( x, y ) the horizontalcoordinates and x in the direction of propagation. Momentum equations, without pressure or gravity read: ( ∂u∂t + u ∂u∂x + v ∂u∂y = ǫ ∇ u ∂v∂t + u ∂v∂x + v ∂v∂y = ǫ ∇ v, (7)with ǫ a viscosity coefficient. Note that mathematical models with a similar formalism have been proposed [3] todescribe the formation of large scale structures in the Universe, assuming a potential flow, namely that a functionΦ( x, y ) exists such that u = − ∂ Φ ∂x and v = − ∂ Φ ∂y . Such a dependence of the velocity field on a potential is possible forthe present model: if the flow is potential at t = 0, equations (7) reduce to the unsteady Bernoulli equation, whichensures that the flow remains potential for all times. In this potential flow case, equations with the viscous termscan still be solved thanks to Florin’s [5] transform (often referred to as Hopf-Cole transform) : S = − ǫ ln | Φ | , whichyields a linear diffusion equation for S ( x, y, t ).For ǫ = 0 equations(7) have the formal solution very similar to equation (3): (cid:26) x − tu ( x, y, t ) = F ( u, v ) y − tv ( x, y, t ) = G ( u, v ) (8)where the functions F, G are presented below. Let M be the general two-by-two matrix of the linear group, withnon-zero determinant. Therefore (˜ u, ˜ v ) = M ( u, v ) and (˜ x, ˜ y ) = M ( x, y ) are solutions of equations of the same formas the left-hand side of Equation (8). This holds true with the same matrix acting on ( u, v ) and ( x, y ). The samegeneral linear map will change in general the functions F ( u, v ) and G ( u, v ) in a rather complicated way. Howeverbecause the map is linear it will not change the degree of a polynomial in u and v . Based on that property, we shallbe able to show, almost without any calculation, that singularities yielding discontinuities of derivatives with respectto x may spread in any direction, the analysis being made by assuming first that the singularity spreads normal to x .The functions F ( u, v ) and G ( u, v ) are defined by the initial conditions. For smooth initial data they are suchthat the mapping from ( u, v ) to ( x, y ), defined by x = F ( u, v ) and y = G ( u, v ) is one to one. Therefore the Jacobian J ( t = 0) = ∂G∂v ∂F∂u − ∂G∂u ∂F∂v never vanishes. The solution of equations (8) becomes singular whenever the time-dependentJacobian J ( t ) of the mapping of ( u, v ) onto ( x, y ) becomes singular. This Jacobian reads : J ( t ) = (cid:18) ∂G∂v + t (cid:19) (cid:18) ∂F∂u + t (cid:19) − ∂G∂u ∂F∂v . (9)Let us assume that the singularity first occurs at t = 0 and x = y = 0, with u = v = 0 as well. The latter restrictionis not as important as one could believe at first : the original equations are Galilean invariant, so that one can alwaysspecify that the velocity field is u = v = 0 at some point of time and space.The condition equivalent to the cancelation in 1D of the second derivative of f ( u ) at u = 0, here, is the propertythat, near t = u = v = 0, the first relevant terms in the TSE of the Jacobian J ( t ) read J ( t ) = at + bu + cv + 2 duv, (10)where a , b , c and d are constant derived from the TSEs of F ( u, v ) and G ( u, v ) near zero. The expansion (10) impliesthat there is no term linear with respect to u and v . Otherwise the Jacobian would vanish for t negative and for smallvalues of u and v , which is not the case based on our assumptions. The two coefficients of the part of the Jacobianlinear with respect to u and v vanish generically at points in the ( x, y ) plane. Of course, besides this cancelation,there are constraints on the sign of a and of the quadratic form bu + cv + 2 duv to make the singularity appear astime increases. For instance, if a is positive, the quadratic form bu + cv + 2 duv must be negative definite, whichrequires c and b negative and d − ab negative.The functions F ( u, v ) and G ( u, v ) have to be expanded in Taylor series at least to third-order to yield the quadraticterm bu + cv + 2 duv in the Jacobian. This depends on two sets of eight coefficients. Fortunately this calculationis greatly simplified by limiting oneself to the quantities at leading order near t = 0. The Jacobian at leading orderdepends on the linear terms in the TSE of F and G near u = v = 0. Let us note that coordinates x and y do nothave to be taken in directions orthogonal to each other, because the Jacobian matrix is reduced to its diagonal formby a general change of coordinates that is not necessarily a rotation, the Jacobian matrix having no reason to be realsymmetric in general. At t = 0 the Jacobian matrix has a zero eigenvalue, that can always be associated with aneigenvector pointing in the x direction. Therefore the coefficients of the TSE of F ( u, v ) and G ( u, v ) have to satisfythe condition that ∂F∂u = 0, ∂F∂v = 0 and ∂G∂u = 0 at u = v = 0. This yields that, near the singularity, one has at leadingorder (cid:26) F ( u, v ) = b ′ v + d ′ u G ( u, v ) = c ′ v, (11)with b ′ , c ′ and d ′ coefficients depending on the initial conditions for the original equation, which are related to theconstant a, b of expansion (10) by the relation a = c ′ , b = 3 c ′ d ′ . Therefore at leading order, the second equation (8)writes y = c ′ v, (12)because one can neglect the term vt as compared to c ′ v (because t is a priori much smaller than the constant c ′ ,assumed non-zero). The first equation (8) becomes: x − b ′ c ′ y − ut = d ′ u . (13)This equation is not yet the one we are aiming for, because it describes a singularity all along the parabola ofCartesian equation x − b ′ c ′ y = 0 at time 0. Actually we omitted a term of order uv in the TSE of F ( u, v ). Let usthus write F ( u, v ) = b ′ v + d ′ u + g ′ uv . The term g ′ uv ultimately yields a term of order y u in the equation for u .Contrary to the term − b ′ c ′ y just considered, this one changes the singularity completely in the sense that, with thisterm, the singularity occurs first at a point in the ( x, y ) plane and later spreads in the y direction. Thus, after someelementary rescaling, equation (6) for the 1D case is changed into:0 = u + x − βy − u ( t − y ) (14)which is the main result of this section.Note, we neglected the term m ′ u v in the TSE of F ( u, v ), although it is of the same order of magnitude, | t | / , asthe one we kept. This is because this term yields a term like u y in the equation for u , which can be eliminated byredefining u as u + m ” y , with m ” = m ′ , a constant. Also, first, note that for actual overturning waves in shallowwater, at least at leading order, the velocity component u is proportional to the surface elevation (see section 4 wheresimple waves are defined for this case). From now on, let us consider u as the surface elevation, until we can furtherdefine the elevation origin. Second, note that in equation (14), the term − βy changes the support of the singularityfrom a straight line to a bent line. Additionally, except for this term, all other terms in equation (14) are of order | t | / near the singularity, with y ∼ | t | / , u ∼ | t | / , x ∼ | t | / . (15)Now, for a given value of t , the extension of the overhanging region of the 3D surface u ( x, y, t ) may be visualized bythe set of points where the tangent plane is vertical; this set of points will be referred to as the apparent contour C v .Differentiating equation (14), we get the relation (3 u + y − t ) du + dx + (2 u − β ) dy = 0. Along C v , the conditionof a vertical tangent, i.e., infinite values of partial derivatives ( ∂u∂x , ∂u∂y ), requires that the coefficient of the du termvanish. Introducing the relation 3 u + y − t = 0 into equation (14), yields the equation for the apparent contour, x = βy ± / ( t − y ) / . (16) - - - xy - - - - xy FIG. 1: Boundary of the overhanging region or apparent contour C v , given by equation (14), with β = − .
5. Leftward figure: just after the beginning of the overturning process at t = 0 .
13; rightward figure : for t = 1 .
4, which is actually beyond therange of validity of the ’generic equation’ .
Equation (16) is plotted in Fig. (1), and shows the boundary of the overhanging domain projected on the horizontalplane, at time t after the inception of the singularity. This boundary spreads along the parabolic curve of Cartesianequation x = βy . The width of the crescent shaped boundary expands as x ∼ | t | / , while its length expands as y ∼ | t | / .These analytical predictions are compared in the next section to results of numerical simulations of wave-breaking,derived from those reported in [6]. C. Numerical simulation of wave-breaking
Euler equations with fully nonlinear free surface conditions have recently been solved numerically in 3D, to simulatethe early stages of solitary wave-breaking induced by changes in topography in shallow water, namely a curved slopingridge [6]. Note that the word ”solitary” refers here to a single, isolated wave, not to a soliton-like wave propagatingwithout changing shape. The flow in the initial solitary wave being irrotational, Kelvin theorem implies that this willbe the case for any positive value of t . Hence, the problem is described by fully nonlinear potential flow equations,which are solved in the model with a high-order boundary element method, combined with an explicit time-integrationmethod, expressed in a mixed Euler-Lagrangian formulation. In fact, earlier comparisons of 2D numerical results withlaboratory experiments have shown that the full potential theory predicts well wave overturning in deep water, as wellas wave shoaling and overturning over slopes [8]. Recent contributions to 3D simulations of nonlinear and overturningwaves are reviewed in [6], where new results for the velocity and acceleration fields during wave overturning are alsocomputed, both on the surface and within the flow. (a)(b) yx z FIG. 2: Numerical results at t = 1 .
378 (based on [6]): (a) 3D wave profile; and (b) apparent contour.
Some of the numerical results reported in [6] are revisited below, and compared to the theoretical predictions ofsection II B.An idealized geometry was selected to simulate 3D breaking-waves on beaches, which is shown in figure 2 of ref.[6].The wave propagates towards the positive x direction, where the water depth z = − h is constant in the first partof the tank, then a sloping ridge starts at x = 5 .
22 with a 1/15 slope in the middle cross-section and a transversemodulation of the form sech ( ky ), with k = 0 .
25 in the present case. (the non-dimensional variables ( x, y ) beingscaled with h and time with p h /g ). On such a bathymetric lens, that causes directional wave focusing, the wavefirst overturns on the axis, at y = 0 and then breaking propagates laterally in both ± y directions. 3D wave profiles areshown in Figs. 5 c, c , and 20, of reference [6], for time values after overturning t = 1 . , .
252 and 1 . t = t ′ − t ′ B in [6]’s notations).In Figure 2a, we reproduce the 3D wave profile at time t = 1 .
378 just before the plunging jet reaches the underlyingsurface of the wave, and show in Figure 2b the apparent contour of the corresponding overhanging region. Figs. 3displays apparent contours projected on the horizontal plane, both for this and an earlier time case, t = 0 .
13. Bothof these apparent contours are in good qualitative agreement with the theoretical predictions of equations (14)-(16),shown in Figure 1 for the same two values of t . Hence, even though the scaling laws were derived above for smallvalues of ( t, u, x, y ) we observe that the ’generic equation’ derived in the previous subsection remains correct beyondthis range of validity. Indeed the apparent contours for the latter time are very similar, see Figs. 2b and 3 (right),while the corresponding values of ( t, u, x, y ) become larger than unity.Let us now test the scaling laws for the longitudinal and transverse expansions in time of the singularity (definedas h ( x, y, t ) being multivalued in the overhanging region), representing the overturning wave. Thus, Figs. 4 and 5show numerical results for the spatial expansion of the apparent contour along x and y , respectively, as a function of x y 17.6 17.8 18 18.2 18.4 18.6−3−2−10123 x y FIG. 3: Numerical results (based on [6]): Projection of the apparent contour at time (left) t = 0 .
13, (right) t = 1 . t t x FIG. 4: Numerical results (+) (based on [6]), compared to analytic prediction, equation (15): Growth of the overhanging regionalong the propagation direction x , as a function of t = t ′ − t ′ B . the time increase from the beginning of the overturning. The two curves display a good fit with the power laws t / and t / , respectively, as predicted by our analysis in section II B. The discontinuity near time t = 0 .
23 in Figure 5 islikely a numerical artifact due to the difficulty of accurately defining the apparent contour along the y direction forsmall values of t . In summary the predictions of the ‘generic solution’(16), proposed above, are well confirmed by thesimulations.While these simple scaling laws have been confirmed, the asymmetry of the actual wave profile, due to the appearanceof a plunging jet, cannot be described by the simple cubic term in the ( F, G ) expansion. In fact, the plunging jet canbe viewed as a late evolution of a Rayleigh-Taylor instability, which occurs when a heavy liquid is placed over a lighterone. This phenomenon was recently studied in [9], where a finger in free fall was also observed in the late (and highlynonlinear) stage of the Rayleigh-Taylor instability. The 2D analysis of [9] assumed a jet falling vertically (along z inour notations) without horizontal velocity (along x ). We will extend this work to estimate further properties of theoverturning jet, by including an initial longitudinal velocity of the jet, u . Thus, in the overturning region, the tip ofthe jet follows a ballistic motion, with z ∼ − gt , and a quasi-constant horizontal velocity (see, e.g., Figure 12a of[6]). Scaling both time and position, as in [9], by factors √ gk and k respectively (where k is the wave number of theimplemented sine perturbation of the interface) and changing their notations y, v into − z, w to avoid confusion withthe above notations, the velocity field in the tip region is approximated by t t y FIG. 5: Numerical results (+) (based on [6]) compared to analytic prediction, equation (15): Growth of the overhanging regionin the direction y transverse to propagation as a function of t = t ′ − t ′ B . (cid:26) u ∼ u w ∼ − t ∼ −√ z. (17)Now, the kinetic equation for the interface α ( x, t ) reads ∂α∂t + u ∂α∂x = −√ α. (18)Defining α ( x, t ) = − t [ t/ − γ ( x, t )], we obtain after linearization ∂γ∂t + u ∂γ∂x = 0 . (19)which can be solved in the form, γ ( x, t ) = θ ( x − u t ). The curvature of the plunging tip is then defined as κ = ∂ α∂t | x = u t .This yields the scaling relationship for the asymptotic evolution of the curvature κ ∼ | t | θ ”0 . (20)[Note that, without the horizontal velocity, one finds κ ∼ | t | θ ”0 .]We now posit that the predictions (17) and (20) of the asymptotic linear evolution of w and κ with time, are alsovalid in the 3D case. Indeed the first one simply reflects the ballistic motion, and the second one remains valid whenthe curvature of the wavefront along the crest is taken into account (this being time independent at leading order, seeequation (16)). Figs. 6 and 7 confirm that the relationships (17) and (20) agree well with results of 3D simulations. III. WAVE-BREAKING IN SHALLOW WATER
The analysis we presented before does not directly apply to actual waves, which are not physically described by thepressure-less fluid equations used so far. However, we will show below that, for surface waves propagating in shallowwater, the equations governing wave motion reduce to the equations of a pressure-less fluid near the singularity.Weakly nonlinear and dispersive inviscid long waves are well described by Boussinesq equations, when
A/h ∼ ( kh ) ≪ A a characteristic wave amplitude, k = 2 π/λ , the wavenumber with λ the wavelength, and h the water depth from mean water level. More recently, extended Boussinesq models (BM) have been proposed thatapply to strongly nonlinear waves, with A/h = O (1), and have been shown to stay valid close to the breaking point[12]. In BMs, the fundamental quantities are the total fluid depth, h ( x, y, t ) (i.e., from free surface to bottom) andthe two Cartesian components of the horizontal fluid velocity, u ( x, y, t ) and v ( x, y, t ). The BM equations of motion t t w t FIG. 6: Numerical result (+) (based on [6]), compared to analytic prediction, second equation (17): Growth of the verticalvelocity at the jet end, at y = 0. . tt κ FIG. 7: Numerical result (+) (based on [6]), compared to analytic prediction, equation (20): Growth of the curvature at thejet end, in the plane y = 0. . express the conservation of mass and of horizontal momentum. Assuming fully nonlinear non-dispersive long waves,BM equations yield the ”nonlinear shallow water wave” equations [11]: ∂h∂t + ∂ ( uh ) ∂x + ∂ ( vh ) ∂y = 0 ∂u∂t + u ∂u∂x + v ∂u∂y + g ∂h∂x = 0 ∂v∂t + u ∂v∂x + v ∂v∂y + g ∂h∂y = 0 , (21)where g is the positive acceleration of gravity. Whenever all functions are independent of y , this set of equationsadmits nontrivial simple 1D wave solutions. We shall see that similar results also exist in 2D. A. One-dimensional case
When the solution only depends on x , we can define h = h + ˜ h ( x, t ) and, if we further assume that h is muchlarger than the free surface perturbation ˜ h ( x, t ), one may neglect[11] terms of order ˜ h ( x, t ) u with respect to terms of0order h u ( x, t ) in the BMs. This yields: ( ∂ ˜ h∂t + u ∂ ˜ h∂x + h ∂u∂x = 0 ∂u∂t + u ∂u∂x + g ∂ ˜ h∂x = 0 . (22)Note that, although ˜ h has been neglected with respect to h , one has kept the term u ∂ ˜ h∂x in the first equation, becausea priori it is not smaller than any other term written explicitly in the second equation (22).Let us now scale the lengths with h and time with p h /g , that amounts to divide the velocity by the long wavecelerity C = √ gh . In terms of scaled variables, the equations for the new functions ˜ h and u thus become: ( ∂ ˜ h∂t + u ∂ ˜ h∂x + ∂u∂x = 0 ∂u∂t + u ∂u∂x + ∂ ˜ h∂x = 0 . (23)If one assumes a simple wave solution of Equations 23, namely that u = ˜ h , the two equations become identical andyield the familiar (so-called) Burgers equation (1), when replacing u (resp. ˜ h ) by ( u + 1) (resp.(˜ h + 1)). This showsthat the singularity in these equations is the same one as found before. This is of course not a new result, because ofthe possibility of solving the full set of equations (21) by Riemann’s method with v = 0 and no dependence on y . B. Waves in two dimensions
Let us now use the full 2D equations (21) and thus define h = h + ˜ h ( x, y, t ), with ˜ h still being small, compared to h . Using the same scaling as above, the equations (21) for the mass and momentum conservation write in terms ofscaled variables: ∂ ˜ h∂t + (cid:16) ∂u∂x + ∂v∂y (cid:17) + u ∂ ˜ h∂x + v ∂ ˜ h∂y = 0 ∂u∂t + u ∂u∂x + v ∂u∂y + ∂ ˜ h∂x = 0 ∂v∂t + u ∂v∂x + v ∂v∂y + ∂ ˜ h∂y = 0 . (24)To check whether these equations approach a pressure-less Bernouilli equation near the wave-breaking singularity,one assumes that the derivative ∂v∂y is non-singular near this singularity, although ∂u∂x diverges (which is true for simplewaves, see below). Therefore one can neglect the former with respect to the latter. Hence, near the singularity, thefirst equation (24) may be replaced by ∂ ˜ h∂t + u ∂ ˜ h∂x + v ∂ ˜ h∂y + ∂u∂x = 0 . (25)Now, if one assumes a simple wave solution, as in the 1D case, namely that u = ˜ h , Equation (25) becomes identicalto the second equation (24). The first and second equations thus coalesce into ∂u∂t + u ∂u∂x + v ∂u∂y + ∂u∂x = 0 . Finallydefining ˜ u = u + 1, the set of equations (24) writes, ( ∂ ˜ u∂t + ˜ u ∂ ˜ u∂x + v ∂ ˜ u∂y = 0 ∂v∂t + ˜ u ∂v∂x + v ∂v∂y + ( ∂ ˜ u∂y − ∂v∂x ) = 0 , (26)which is identical to the pressure-less fluid system (7), because the expression ( ∂ ˜ u∂y − ∂v∂x ) vanishes for irrotationalfluids.This shows that, at least at leading order, wave-breaking scaling laws are identical for a pressure-less fluid and forthe nonlinear shallow water wave equations.This brings an interesting point: it seems that the widening of the wave-breaking domain in the y direction, likethe square root of time, is ’universal’ for non-linear non-dispersive waves. Other kinds of waves are known to break,like ocean waves on deep water (with, practically infinite depth). Does their spreading in the y direction follow thesame square root law ? The answer to this question depends in particular on how precisely one can identify the timefor the onset of the breaking process. This deserves further investigations.1 IV. DISCUSSIONA. Regularisation
The regularization of the singularity appearing in solutions of the inviscid Burgers equation has been extensivelystudied. One such regularization scheme relies on an assumed balance of nonlinear and dispersive effects over a non-vanishing fluid depth. This yields the well-known Korteweg-de Vries (KdV) equation in 1D. This equation has beenvery thoroughly studied and can be shown to be a form of the standard Boussinesq equations when assuming wavepropagation occurs in a forward direction only [11]. Hence, as for the standard weakly nonlinear and dispersive BMmentioned above, KdV equation requires that the two physically unrelated processes of nonlinearity and dispersionbe of the same order. Hence, the range of applicability of this equation to water waves is quite narrow.A second way to avoid the singularity consists in adding a small dissipative term to the inviscid Burgers equation(1) to make it ressemble the Navier-Stokes equation. This yields the so-called Burgers equation, ∂u∂t + u ∂u∂x = ǫ ∂ u∂x , (27)with the coefficient ǫ being small and positive. With a non-zero ǫ value, the solution of equation (27) is not multi-valued anymore (the equation is regularized) and one has to replace in the above analysis the multi-valued U by afunction with a jump at ζ = 0, which represents a shock wave for this model. By symmetry, the shock discontinuityis at ζ = 0. There, U jumps from +1, for slightly negative values of ζ , to ( − ζ . In asmall interval near ζ = 0, of width of order ( ǫt ) / , the solution goes continuously from +1 to ( −
1) across the shockfront. Moreover, just at the time of the discontinuity, there is also a time interval where the added diffusion term isnowhere negligible. This time interval is such that the thickness of the self-similar region in x is of the same order asthe typical length scale associated with diffusion. The former decreases to zero as t tends to zero, like | t | / , and thelatter like ( ǫt ) / . The cross-over time is such that the two length scales are of the same order of magnitude, whichhappens for t ∼ ǫ . For time scales of this order or shorter, the self-similar solution is not valid anymore and anotherlimit has to be taken. Real wave-breaking, however, cannot be described by this regularization, because it implies anoverturning of the surface that is not represented by the viscous-like term on the right-hand side of equation (27). B. Geometrical meaning of the singularity or where is the true singularity?
Before proceeding further, it is worthwhile reconsidering the question of the singularity in light of a geometricalapproach. This will yield a more precise definition of what we mean by wave-breaking. Let us get back to the ’generic’solution (6) of the pressure-less fluid equation (1). In the simple wave case, this solution takes a simple geometricalmeaning in which the two quantities u (velocity component along the propagation direction) and h (elevation of thefluid surface) are proportional to each other, up to irrelevant multiplicative constants. Then the 1D solution of,0 = h + x − ht (28)can be seen as describing the occurrence of an apparent contour of a curve drawn on a plane, as it is continuouslydeformed, the deformation parameter being t . Such an apparent contour is defined (see section II) as the set of pointswhere the tangent is parallel to the direction of observation, and the time origin t = 0 is the instant where all pointscollapse into a single one. It is worthwhile noting that this definition depends upon the direction of observation.Indeed in 1D, at t = 0, the local Cartesian equation of the curve is h + x = 0, h being the abscissa and x theordinate. As time changes, the curve is deformed and generically takes a non zero slope at the origin. The slopeincreases linearly with time, while the two extrema of the curve x ( h, t ) move apart from each other, as ∆ x ∼ t / .Therefore the equation (28) describes the local shape of a line, when it acquires a contour at t = 0. This contouris fundamentally linked to the direction of observation and would vanish if this direction were tilted in the directionparallel to the slope, at the inflection point. While it seems natural to choose the vertical axis as the observationdirection in real life, we must note that the description of wave-breaking as the appearance of a contour is ambiguous.To get round this difficulty let us note that wave-breaking could be defined either by the occurrence of overturning(i.e., an overhanging region on the surface), or by the impact and intersection of the plunging jet with the surface un-derneath. Casual observation indicates that the two phenomena are very strongly linked in the sense that overturningusually leads first to plunging and then to impact with the underlying surface. However things are not generally sosimple, especially because the equations of motion are reversible (viscosity is negligible and thus neglected for thistypically large Reynolds number flow). Hence, under time-reversal symmetry, a plunging jet would first retract ontoitself and later return to a disappearing apparent contour.2Such a reversible event, namely the possibility of reversible overturning (i.e., the overhanging region spontaneouslyretracting instead of becoming a plunging jet), could possibly occur. This could happen, for instance, for a wavepropagating over a decreasing depth, that would initiate overturning, then entering a region of rapidly increasingdepth, which could make the overhanging region gradually disappear, provided its development has not proceededtoo far in time. Another case could be that of a surface pressure (i.e, caused by wind) gradually increasing and‘pushing’ against the nascent overturning region. Therefore an upward moving retracting jet is certainly among thesolutions of the full dynamical problem and would lead quite naturally to define wave-breaking as the occurrence ofself-penetration of the surface, a definition free of arbitrariness because, contrary to anything related to the apparentcontour, it does not rely on an arbitrary choice of the direction of observation.Seen from another perspective, wave-breaking now defined by the self-penetration of the surface by the plungingjet follows, but not always, from overturning of the surface. This raises also the very interesting issue of describingthe post-penetration flow and surface intersection behavior (discussed below). Indeed regularization by capillaryphenomena, viscosity, etc. . . may make this a well posed problem. Nevertheless we posit that the problem of impactbetween the plunging jet and the surface underneath it can be formulated in almost the same manner as the problemof propagation of nonlinear waves. This is discussed in the next subsection. C. Post collision dynamics
Classically the dynamics of nonlinear waves in the inviscid limit can be formulated as a system of two coupledfirst-order (in time) equations for the value of the velocity potential on the free surface and for the position of anarbitrary point on this surface. Let us briefly summarize the bases for this approach, in the usual case withoutcollision of fluid interfaces. Knowing the value of the velocity potential on the free surface and, given the Neumanncondition on the bottom, one can in principle solve the Dirichlet-like Boundary Integral Equation (BIE) problem forthe velocity potential everywhere [6, 13]. This yields in particular the velocity normal to the free surface that is usedto advance this surface in time. The value of the velocity potential on the surface is then similarly advanced by timeintegrating the Bernoulli equation at every free surface point. The major advantage of this BIE method is that, innumerical studies, one can transform the problem of solving Laplace or Poisson’s equation over the entire domain intoan integral equation on the surface, solvable by analytic function theory in 2D and by Green’s function method inboth 2D and 3D [8, 13], short-cutting the computationally complex problem of drawing a smooth curve or a surfaceembedded in the geometrical space.To summarize, one starts from the unsteady Bernoulli equation for the velocity potential Φ( r , t ): ∂ Φ ∂t + 12 ( ∇ Φ) + gz + pρ = 0 (29)with z the height, ρ the mass density of the fluid and p the pressure. This equation is expressed on the free surfaceof the fluid, located at z = h ( r H , t ), where r H is the set of horizontal coordinates. Note that if there is overturning,the function h ( r H , t ) is not single-valued everywhere. Because the free surface is in contact with air, the atmosphericpressure is applied everywhere on the surface, whose variation can be neglected. Therefore, we can set p = p a = 0 at z = h and equation (29) becomes: ∂ Φ ∂t = −
12 ( ∇ Φ) − gz | z = h ( r H ,t ) . (30)The other conditions to be satisfied by Φ are mass conservation, which becomes Laplace’s equation: ∇ Φ = 0 , (31)and the Neumann (no-flow) condition at the bottom (assumed for simplicity at z = 0), ∂ Φ ∂z = 0 | z =0 . (32)Given the value of Φ on the free surface and the Neumann condition at z = 0, one can solve Laplace’s equation, onceit is transformed into a linear BIE. Concretely this means that the normal derivative of Φ on the free surface can befound by solving an integral equation on the boundary, the remaining components of the gradient of Φ, there, beingknown from the data. This allows us to compute the left-hand side of equation (30), from its right-hand side, andthus to advance the value of Φ on the free surface, in time. The evolution of the free surface elevation is given by3the kinematic condition that it be a material surface, thus advected by the local fluid velocity ( w is the vertical fluidvelocity): ∂h∂t + ∂ Φ ∂x ∂h∂x + ∂ Φ ∂y ∂h∂y = w | z = h ( r H ,t ) . (33)Equations (30) and (33), assuming that the velocity components are computed via the solution of Laplace’s equation(31) (by a BIE), define a closed problem for the coupled evolution of the velocity potential and the free surface. Allof this, however, assumes that the free surface remains simply connected. The explicit form of the BIE problem is asfollows. It uses the 3D free space Green’s function defined as [10], G ( x , x ) = 14 πr , (34)with the normal derivative ∂G ( x , x ) ∂n = − π r · n r , (35)and r being the distance from point x = ( x, y, z ) to the reference point x = ( x , y , z ), both being on the boundary,and n representing the outward normal unit vector to the boundary at point x . Green’s second identity transformsLaplace’s equation into the BIE, α ( x )Φ( x ) = Z Γ( x ) (cid:20) G ( x , x ) ∂ Φ( x ) ∂n − Φ( x ) ∂G ( x , x ) ∂n (cid:21) dΓ , (36)in which α ( x ) = θ π , with θ the exterior solid angle made by the boundary at point x (i.e. 2 π for a smoothboundary). Moreover dΓ is the element of area on the surface at position x . The accurate numerical solution of BIE(36), with the bottom boundary condition (32), together with the time integration of Equations (30) and (33), allowssimulating 3D overturning waves such as depicted in Figure 2a [6, 13], up to the time the plunging jet is about toimpact the free surface underneath it.
1. Generalisation of the Integral equation method
In the case of a plunging jet impacting the underlying surface, the domain becomes multi-connected and the aboveBIE cannot be used as such anymore. We discuss a way to extend this method to such cases. According to theKelvin-Helmholtz theorem, for an inviscid fluid, no vorticity can be created outside of the domain boundary. In thepresent case, this applies to the boundary between the free falling overturning jet and the fluid underneath. Thus,assuming the jet does not immediately mixes with the impacted fluid (and this is in fact supported by some laboratoryexperiments for early times after impact) a vortex sheet can be specified on the surface separating the two fluids [14],to simulate the tangential velocity discontinuity that appears, such that the plunging jet may continue to exist as aregion of potential flow beyond the time of intersection. Thus the jet has two parts, an upper one in contact withambient air (where p = 0) , and a lower part inside the underlying fluid, which is separated from the surroundingfluid by high pressure air (see Figure 8). Let us now derive a closed set of dynamical equations for quantities definedon surfaces, as in the previous case prior to intersection. For that purpose, we introduce three boundaries, labeled B , B and B , respectively (see Figure 8). B is the bottom at z = 0, B is the free surface between the air andthe fluid, where the pressure is zero, and B is the inner surface separating the fluid coming from the jet from theunderlying fluid.Because of the existence on B of the vortex sheet, i.e., concentrated vorticity, one cannot take as dynamicalfunction the value of the velocity potential there, since it jumps across this interface. Therefore we shall take ∂ Φ( x ) ∂n instead of Φ( x ), as dynamical variable on B where the normal velocity and the normal acceleration have both to becontinuous, to make the two sides of the boundary moving together. The BIE for Φ can be solved from the knowledgeof ∂ Φ ∂n (Neumann condition) on B , B , and Φ on B (Dirichlet condition). This yields the set of boundary conditionsfor solving the Laplace equation (31) inside the fluid. To advance in time and derive the dynamical equation for ∂ Φ ∂n on B , let us define q = pρ + | ∇ Φ | that we shall call the dynamical pressure. By taking the Laplacian of theBernoulli equation (29), and because Φ and z are harmonic functions, one finds:∆ q = 0 . (37)4 FIG. 8: Scheme of the plunging jet.
We are now ready to solve equation (37), using only two boundary conditions for q , q = | ∇ Φ | on the freesurface B where p = 0, plus the condition ∂ ( q − |∇ Φ | ) ∂z + g = 0 on the bottom B . Indeed the portion B does notenter in the contour of the BIE for q , because both q and ∂q∂n are continuous on B for the following reasons.Defining the velocity as −→ u = −−→∇ Φ , the Euler fluid equation can be written as, ∂ −→ u∂t + −→∇ q + ge z = , (38)where e z is the unit vector in the vertical direction. The continuity of q across B follows from equation (38) sinceotherwise an infinite acceleration would take place. Moreover the continuity of ∂q∂n is ensured because of the continuityof the normal acceleration across B . Therefore the Laplace equation (37) for q can be solved by using the boundaryconditions on B and B only.This means that, like the problem of evolution of the free surface in the simply connected case, the post-collisiondynamical equations (for Φ and q ) can be reduced to linear integral equations posed on the boundaries ( B , , for Φ,and B , for q ).Because of the merging of the plunging jet with the fluid underneath, vorticity is generated in a flow that is otherwisepotential. This vorticity is roughly horizontal in the span-wise direction of the waves. If those waves are more or lessparallel to each other (i.e., long crested) , this could be a coherent source of vorticity near the surface that wouldbecome on a large scale a shear layer driving the fluid opposite to the direction of the wave propagation. Anotherdifficulty in this problem is the possible occurrence of regions of negative pressure where nucleation of bubbles shouldtake place (actually of pressure less than the value of saturating pressure for the given temperature condition). Weset aside this possibility, although it is obvious from observations that air bubbles are generated in white caps andeven more so by plunging breakers.
2. Hertz-fluid contact problem
A possible explanation for this generation of many bubbles lies in a simple scaling law for the Hertz-fluid contactproblem. Let us assume that the plunging jet has a locally parabolic tip hitting the fluid underneath, for simplicity,assumed to be flat. Therefore we have a typical length scale, the radius of curvature of the parabolic tip R , anda velocity scale, the velocity w of the falling jet when it hits the flat surface. Because Bernoulli equation has nointrinsic scale, every parameter, like for instance the typical value of the velocity in the collision region (defined a bitmore accurately below), is like u ( x, t ) = w u sc ( x/R, wt/R ), where the vector field u sc ( x/R, wt/R ) is dimensionless,universal, and defined by numerical functions only, t = 0 being defined as the time of the first contact. One may tryto find the behavior of the various quantities involved just after the collision by borrowing some of the ideas of thefamous Hertz contact problem [15] between solids. As in Hertz contact, one assumes that the perturbation comesfrom the flattening of the tip, wherever it would be inside the fluid underneath if this were possible. The tip wouldhave penetrated by a length wt inside this fluid, and so spread over a length of order λ ∼ ( Rwt ) / along this surface.Because this brings a perturbation of order 1 to the boundary of the fluid inside the jet and because Laplace’s equationhas no typical length scale (this argument is directly borrowed from Hertz’s contact problem, where Hertz used this5property of the equations of elasticity), the region perturbed inside the jet is a volume of size λ in all directions. Thevertical momentum of the falling jet in this perturbed region is therefore of order P z ∼ ρwλ . (39)The change with time of this momentum dP z dt = pλ is due to the pressure p near the boundary between the fallingjet and the fluid underneath over an area λ . This yields an order of magnitude for the pressure in the collisiondomain p ∼ ρw (cid:18) Rwt (cid:19) / , (40)when using the approximation dP z dt ∼ P z t . The relation (40) shows that this pressure, felt across the boundary B ,diverges at the instant of contact where it is regularized probably by capillary phenomena and/or compressibilityeffects in the fluid (in the compressible case, at short time, the divergence of pressure is replaced by a very large peakof pressure of order ρwc s , with c s the speed of sound in the fluid). This result could explain in part the formation ofbubbles and foam in the real world. V. FINAL REMARKS
We presented a rather simplistic, but physically meaningful, approach to wave-breaking, yielding scaling laws thatwere shown to be in good agreement with numerical results for wave overturning in 3D, based on a higher-ordersolution of full dynamic equations. This work is an example of complementarity between advanced numerical andanalytical methods, each providing validation and physical meaning to the other. Analytical methods also make itpossible to generalize nonlinear properties observed in numerical results, to a whole class of equations and problems.In particular, and on a more general point of view, it seems to us that many important concepts and ideas pertainingto nonlinear science remain to be exploited, in many areas where they can guide the intuition or better yield specificpredictions that can then be compared to numerical results. In the specific application of nonlinear fluid mechanicspresented here, we have encountered various more generic subproblems, all with a definitely very strong nonlinearflavor, like the singular behavior of solutions of partial differential equations, or what we referred to as the Hertz-fluidcontact problem. In closing, this indicates that free boundary problems remain a rich and perhaps inexhaustiblesource of inspiration in the field of nonlinear science.
Acknowledgments
Paul Clavin, Christophe Josserand, Laurent Di Menza and Fr´ed´eric Dias are gratefully acknowledged for stimulatingdiscussions. Philippe Guyenne acknowledges support from the University of Delaware Research Foundation and theUS National Science Foundation, under grant DMS-0625931. St´ephan Grilli acknowledges support from the U.S.ONR Coastal Geosciences Division (code 321CG) Grant No. N000140510068. [1] See for example the review paper by Dias F and Kharif C 1999 Nonlinear gravity and capillary-gravity waves
Annu. Rev.Fluid Mech. Journal de l’Ecole Polytechnique
14 ´eme cahier Astrofisica Astrophysics Catastrophe theory (Springer Verlag, Berlin). See in particular Chapter 8 on Caustics and wavefronts.[5] Florin V A 1948
Izvestia Ak.Nauk SSSR Otd.Tekhn.Nauk J. Fluid Mech.
Proc. R. Soc. A
Compte-Rendus de la rencontre du Non-lin´eaire 2008,Paris (Non-linaire Publications, Orsay France). [8] Dommermuth D G, Yue D K P, Lin W M, Rapp R J, Chan E S and Melville W K 1988 Deep-Water Plunging Breakers :a Comparison between Potential Theory and Experiments J. Fluid Mech.
J.Fluid Mech.
J.Waterway Port Coastal Ocean J. Waterway Port Coastal Ocean Engng.
Phys. Rev. Lett. The Boundary Element Method for Engineers (Wiley, New York) .[11] Mei C C 1983
The Applied Dynamics of Ocean Surface Waves ( World Scientific, Singapore).[12] Wei J, Kirby J T, Grilli S T, and Subramanya R A 1995 A Fully Nonlinear Boussinesq Model for Surface Waves. Part 1.Highly Nonlinear Unsteady Waves
J. Fluid Mech.
J. Wtrwy, Port, Coast, and Oc. Engrg. Advancesin Coastal Modeling ( V. C. Lakhan (ed) , Elsevier Oceanography Series ).[13] Grilli S T, Guyenne P, and Dias F 2001 A fully nonlinear model for three-dimensional overturning waves over arbitrarybottom