Next Article in Journal
Operational Modeling of North Aegean Oil Spills Forced by Real-Time Met-Ocean Forecasts
Next Article in Special Issue
Ship’s Wake on a Finite Water Depth in the Presence of a Shear Flow
Previous Article in Journal
Wave Forces on a Partially Reflecting Wall by Oblique Bragg Scattering with Porous Breakwaters over Uneven Bottoms
Previous Article in Special Issue
Wave Energy Dissipation of Spilling and Plunging Breaking Waves in Spectral Models
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A 2D Model for 3D Periodic Deep-Water Waves

by
Dmitry Chalikov
1,2
1
Shirshov Institute of Oceanology RAS, 117997 Moscow, Russia
2
Department of Infrastructure Engineering, Melbourne School of Engineering, University of Melbourne, Victoria 3010, Australia
J. Mar. Sci. Eng. 2022, 10(3), 410; https://doi.org/10.3390/jmse10030410
Submission received: 13 February 2022 / Revised: 4 March 2022 / Accepted: 7 March 2022 / Published: 11 March 2022
(This article belongs to the Special Issue Nonlinear Wave Dynamics and Wake Structure)

Abstract

:
The paper is devoted to further development of an accelerated method for simulation of the two-dimensional surface waves at infinite depth with the use of a two-dimensional model derived with simplifications of the three-dimensional equations for potential periodic deep-water waves. A 3D full wave model (FWM) is based on a numerical solution of a 3D Poisson equation written in the surface-fitted coordinates for a nonlinear component of the velocity potential. For sufficient vertical resolution used for the Poisson equation, the 3D model provides very high accuracy. The simplified model is based on the 2D Poisson equation written for a free surface. This exact equation contains both the first and second derivatives of the velocity potential, i.e., it is unclosed. The analysis of the accurate solutions for the 3D velocity potential obtained with the 3D model shows that those variables are linearly connected to each other. This property allows us to obtain a 2D equation for the first derivative of the velocity potential (i.e., vertical velocity on surface), which gives the closed 2D formulation for a 3D problem of two-dimensional waves. The previously developed scheme was not universal since the parameters of the closure scheme had to be adjusted to the specific setting. The current paper offers a new formulation of a closing scheme based on the integral parameters of wave field. The method of closing the equations, as well as the numerical parameters, were chosen on the basis of the multiple numerical experiments with the full nonlinear wave model (FWM) and selection of a suitable closing scheme. That is why the given model can be called Heuristic Wave Model (HWM). The connection between the first and second variables is not precise; hence, the method as a whole cannot be exact. However, the derived 2D model is able to reproduce different statistical characteristics of the 2D wave field with good accuracy. The main advantage of the model developed is its high performance exceeding that of 3D model by about two decimal orders.

1. Introduction

A brief review of different 3D numerical methods developed for investigation of wave processes is given in [1]. Among the most advanced approaches, the following methods should be mentioned: (1) the Boundary Integral Method [2,3,4]; (2) a grid model designed at Technical University of Denmark (see [5]); (3) a High-Order Spectral (HOS) model developed at Ecole Centrale Nantes, LHEEA Laboratory [6]. A 3D model for direct simulation of periodic waves combining Fourier transform method with the finite-difference approximation in vertical for Poisson’s equation was suggested in [7] (see also [8]). The last model (we call it “Full Wave Model”, FWM) was used for simulation of wave field development under the action of the wind [9]. In the current paper, FWM for infinite depth is used as a basic model for building up a simplified model that allows for obtaining similar results at a much lower cost.
The main drawback of the 3D phase-resolving models is their low performance because all of them, one way or other, reproduce the vertical structure of wave field on the basis of a 3D equation for the velocity potential. For example, a model [1] with a resolution of 512 × 513 modes ( 2048 × 1024 knots) and just 10 vertical levels used for the simulation of wave field development over 1,200,000 time steps was running in one medium speed processor for about 45 full days. However, especially at the stage of development, the model should be run hundreds of times for validation of schemes and tuning of the parameters. Instead, just several relatively short runs before a long numerical experiment mentioned above were conducted. This is obviously not enough.
It is well known that most of the computer time for a 3D model is spent calculating the 3D structure of the velocity potential. It is necessary to emphasize that the 3D solution is used for calculation of the vertical velocity w on a surface only. Unlike the 2D conformal model, a simple method of calculation of w for a 3D case is not known. In the current paper, a new method of approach to the problem is suggested. The method is based on separation of the velocity potential φ into the linear φ ¯ and nonlinear φ ˜ components and on consideration of Poisson’s equation for the nonlinear component on a surface. This exact equation contains both the first φ ˜ ζ = w ˜ and second φ ˜ ζ ζ = w ˜ ζ derivatives, which means that the formulation of the entire problem is not closed. However, the multiple numerical calculations with an accurate 3D model [9] showed that the vertical profiles of a nonlinear component of the velocity potential φ ˜ have quite a universal simple structure, which allows the suggestion that w ˜ and w ˜ ζ can be closely connected. It was shown in (10) that dependence of w ˜ on w ˜ ζ is very close to the linear one, but the coefficients in such connection depend on external integral parameters of the wave field. In the current paper, a new connection between w ˜ and w ˜ ζ with the coefficients depending on nondimensional parameters is tested. This suggestion turned out to be reliable enough to formulate the entire model that can realistically simulate wave dynamics in a broad range of external parameters.
Note that parameterization was made based on numerical experiments with 3D. FWM for infinite depth; hence, the 2-D model is suitable only for the case of infinite depth.

2. Derivation of a Simplified 2D Wave Model (HWM) for 2D Waves

The equations are written in the non-stationary surface-following non-orthogonal coordinate system:
ξ = x , ϑ = y , ζ = z η ξ , ϑ , τ , τ = t
where η x , y , t   = η ξ , ϑ , τ is a moving periodic wave surface defined by Fourier series
η ξ , ϑ , τ   = M x < k < M x M y < l < M y h k , l τ Θ k , l
where k and l are the components of the vector wave number k , h k , l τ are Fourier amplitudes for elevations η ξ , ϑ , τ , M x and M y are the numbers of modes in the directions ξ and ϑ , respectively, whereas Θ k , l are the Fourier basic functions represented as a matrix:
Θ k l = cos ( k ξ + l ϑ ) 1 k M x , M y l < M y cos ( l ϑ ) k = 0 , 0 l M y sin ( l ϑ ) k = 0 , M y l 1 sin ( k ξ + l ϑ ) M x k 1 , M y l < M y
The surface conditions for potential waves in a system of coordinates (1) at ζ = 0 take the following form:
η τ = η ξ φ ξ η ϑ φ ϑ + 1 + η ξ 2 + η ϑ 2 φ ς
φ τ = 1 2 φ ξ 2 + φ ϑ 2 1 + η ξ 2 + η ϑ 2 φ ζ 2 η p
where φ is the velocity potential. Laplace equation for φ at ζ 0 turns into the Poisson equation
φ ξ ξ + φ ϑ ϑ + φ ζ ζ = ϒ φ ,
where ϒ is the operator:
ϒ ( ) = 2 η ξ ( ) ξ ζ + 2 η ϑ ( ) ϑ ζ + η ξ ξ + η ϑ ϑ ( ) ζ η ξ 2 + η ϑ 2 ( ) ζ ζ .
Strictly speaking, Equations (4)–(6) should be solved jointly, but since the Fourier amplitudes h k , l τ   = h k , l 2 + h k , l 2 1 / 2 for elevation and φ k , l τ   = φ k , l 2 + φ k , l 2 1 / 2 for the surface velocity potential change slowly, some sort of time-splitting scheme is used: Equations (4) and (5) are used as evolutionary equations for η and φ , whereas Equations (6) and (7) are used to satisfy the continuity conditions. It is accepted that in the process of the solution of Equation (7) the coefficients depending on η are constants. Introducing the correction of η and φ s with Equations (4) and (5) in the course of solution does not provide any advantages.
Equations (4)–(7) are written in a non-dimensional form by using the following scales: length L, where 2πL is a dimensional period in a horizontal direction; time L / g 1 / 2 (g is the acceleration of gravity) and the velocity potential L 3 / 2 g 1 / 2 . The pressure is normalized by water density so that the pressure scale is Lg. Equations (4)–(6) are self-similar to transformation with respect to L. The dimensional size of the domain is 2 π L × 2 π L M y / M x . All of the results presented in this paper are nondimensional.
It is suggested in [7] that it is convenient to represent the velocity potential as a sum of two components: an analytical (‘linear’) component φ ¯ and an arbitrary (‘non-linear’) component φ ˜ :
φ = φ ¯ + φ ˜ .
The analytical component φ ¯ satisfies Laplace equation:
φ ¯ ξ ξ + φ ¯ ϑ ϑ + φ ¯ ζ ζ = 0 ,
with the known solution:
φ ¯ ( ξ , ϑ , ζ , τ ) = k , l φ k , l ( τ ) exp k ζ Θ k , l ,
where k =   k 2 + l 2 1 / 2 , φ k , l are Fourier coefficients for the surface potential φ ¯ at ζ = 0 . The 3D solution of (9) satisfies the following boundary conditions for the surface and at infinite depth:
ς = 0 : φ ¯ 0 = φ ς : φ ¯ ζ 0
The presentation (8) is not used for solution of evolutionary Equations (4) and (5), because it does not provide any noticeable improvements of accuracy and speed of the calculations.
The nonlinear component satisfies an equation:
φ ˜ ξ ξ + φ ˜ ϑ ϑ + φ ˜ ζ ζ = ϒ φ ,
Note that the full potential φ on the surface is used as a boundary condition for a linear component, so, a nonlinear component φ ˜ on the surface is equal to 0. Thus, Equation (12) is solved with the boundary conditions:
ς = 0 : φ ˜ = 0 ς : φ ˜ ζ 0 .
Hence, on the surface ζ = 0 Equation (12) takes the form:
φ ˜ ζ ζ = ϒ φ ˜ + ϒ φ ¯
The derivatives of a linear component φ ¯ in (7) are calculated analytically. A method of solution combines a 2D Fourier transform method in the ‘horizontal surfaces’ and a second-order finite-difference approximation on the stretched staggered vertical grid providing high accuracy in the vicinity of the surface. An Equation (12) in FWM is solved as Poisson equations with iterations over the right-hand side by the TDMA method [10] with the prescribed accuracy. The number of iterations depends on a maximum of the local surface steepness, and for high steepness it does not exceed three.
A 3D model described above is potentially exact. Its only drawback is a low speed of performance. Below is a scheme that at a cost of minor loss of accuracy provides a significantly higher speed of calculation.
Since φ ˜ 0 = 0 , an Equation (12) for the velocity potential at ζ = 0 takes the form:
w ˜ ζ = 2 ( η ξ w ξ + 2 η ϑ w ϑ ) + Δ η w s w ζ .
Here, the symbols Δ = ξ ξ + ϑ ϑ and s = η ξ 2 + η ϑ 2 are used. It is remarkable that an Equation (15) is exact. From Equation (14) it follows that
w ˜ ζ = 2 η ξ w ¯ ξ + w ˜ ξ + 2 η ϑ w ¯ ϑ + w ˜ ϑ + Δ η w ¯ + w ˜ s w ¯ ζ + w ˜ ζ ,
hence, Equation (15) can be represented as follows:
1 + s w ˜ ζ = 2 η ξ w ˜ ξ + η ϑ w ˜ ϑ + Δ η w ˜ + r ¯ ,
where the term r ¯ depends only on a linear component φ ¯ :
r ¯ = 2 η ξ w ¯ ξ + η ϑ w ¯ ϑ + Δ η w ¯ s w ¯ ζ
which is calculated analytically using Fourier presentation (10) by the formulas:
w ¯ ( ξ , ϑ , ζ ) = k , l k φ k , l exp k ζ Θ k , l ,
w ¯ ζ ( ξ , ϑ , ζ ) = k , l k 2 φ k , l exp k ζ Θ k , l
w ¯ ξ ( ξ , ϑ , ζ ) = k , l k k φ k , l exp k ζ Θ k , l ,
w ¯ ϑ ( ξ , ϑ , ζ ) = k , l l k φ k , l exp k ζ Θ k , l
The presence of w ˜ ζ in Equation (17) makes the system of equations unclosed. Previously, the author regarded a surface condition (17) as a means of validation of the numerical scheme for the 3D Full Wave Model (FWM) based on a numerical solution of an equation for the velocity potential (see example of such calculation in [8]). This validation is not trivial since w ζ is not a direct product of the Poisson equation solution. The finite-difference derivatives w ˜ = φ ˜ ζ and w ˜ ζ = φ ˜ ζ ζ after solution of the 3D equation for a nonlinear component of the velocity potential were calculated as follows:
w ˜ = φ ˜ ζ ζ = 0 = D 1 ζ 2 2 φ ˜ 1 ζ 1 2 φ ˜ 2 ,
w ˜ ζ = φ ˜ ζ ζ ζ = 0 = 2 D 1 ζ 2 φ ˜ 1 + ζ 1 φ ˜ 2 .
Here, D = ζ 1 ζ 2 2 ζ 2 ζ 1 2 , ζ i is depth where φ i is located, ζ 0 = 0 , the thickness of layer Δ i = ζ i ζ i + 1 grows exponentially with depth ζ : Δ i + 1 = γ Δ i , whereas the value of Δ 0 is calculated from the condition i = 0 i = L 1 Δ i = H where L is the number of levels, and H = 4 π / k p is depth of the domain ( k p is a wave number of a mode with the largest amplitude).
Since a relation (17) formally follows from the Poisson equation, the agreement between these values characterizes accuracy of the simplified scheme. The agreement cannot be very precise since the variables w ˜ ζ and w ˜ are not a direct product of the solution, and they are calculated by the approximations (23) and (24) following the solution of the Poisson equation. Moreover, the approximations for solution of the Poisson equation and calculation (23) and (24) are different. However, the calculations of w ˜ ζ on the basis of Poisson Equation (12) with the approximations (23) and (24) and calculations by a surface condition (17) showed a very good agreement (see Figure 1 in [9]), which prompted us to carry out further analysis of Equation (18).
It was found that variables w ˜ and w ˜ ζ turned out to be well connected to each other. Firstly, the simplest linear hypothesis was tested:
w ˜ = A 0 + A 1 w ˜ ζ
where A 0 and A 1 are constants. This dependence was considered in [11] with the use of the 3D FWM model (described in [1,7,8,9]) that calculates the 3D structure of the velocity potential as well as 2D fields of w ˜ and w ˜ ζ . The dependence between w ˜ and w ˜ ζ can be studied in Fourier space as well as in a physical space. The results obtained in Fourier space turned out to be far less certain than in a physical space, i.e., the dependence between the Fourier coefficients for w ˜ and w ˜ ζ showed considerable scatter. In our opinion, it happens because the medium and high wave numbers of Fourier amplitudes are generally unstable, whereas the verticals structure of a physical layer involved in orbital movement is universal because it is composed of a large number of modes.
The calculations performed in [11] were repeated here at a higher grid resolution 2048 × 1024 and a number of vertical levels equal to 50. The dispersion of surface elevation σ , as well as the integral steepness s, variate in the ranges 0.033 < σ < 0.066 and 0.02 < s < 0.20 , respectively. The fields included 1.7 · 10 8 pairs of values w ˜ and w ˜ ζ . The statistical connection between those variables is shown in Figure 1. It is important to underline that dependence (25) is valid for the entire wave field, i.e., the parameters A 0 and A 1 are close to constants. Otherwise, the whole approach based on a hypothesis (25) might be incorrect.
The data used for analysis included the 3D velocity potential fields φ ˜ ξ , ϑ , ζ , τ and the values of w ˜ and w ˜ ζ calculated by (23); the values of w ˜ ζ calculated by a relation (17) and some other characteristics (see below). Since a high value of γ = 1.25 for generation of vertical grid was used, the accuracy of calculation of the derivatives (23) and (24) was very high. The comparison with the analytical values of the derivatives shows that the errors of approximation (23) and (24) do not exceed 10 10 .
The surface values of w ˜ and w ˜ ζ are defined by an asymptotic behavior of a nonlinear component of the velocity potential φ ˜ at ζ 0 . Hence, the connection between them can be found in the close vicinity of ζ = 0 .
A dimensional form of (25) is
W = L ( A 0 + A 1 ) W ζ ,
where W and W ζ are dimensional analogues of w ˜ and w ˜ ζ , whereas L is the external scale of the problem. Formally, for any specific wave field which is more or less horizontally homogeneous, the dependence (25) is correct at appropriate choice of coefficients. However, the problem is that for different wave fields (for example, for different stages of wave field development) the variability of those coefficients can be significant. It happens because the scaling with external scale L does not reflect the internal properties of wave field. To provide a more universal approximation, it is necessary to introduce the scaling in terms of integral characteristics of a given wave field with no use of L .
For evaluation of connection between w ˜ and w ¯ ζ , ten thousand short-term numerical experiments with FWM were performed. The Fourier resolution is equal to 257 × 257 modes; the number of levels for solution of Poisson equation equals 30, whereas a stretching coefficient γ = 1.25 . The initial conditions for elevation were calculated on the basis of JONSWAP spectrum [12] for the inverse wave age U / c p = 1 ( U is wind velocity, c p is the phase velocity in peak of spectrum). The calculations were carried out for different values of peak wave numbers in the range 10 k p 100 . The data collected are characterized by the nondimensional dispersion of elevation
σ = η η ¯ 2 ¯ 1 / 2 ,
in the range 0.003 < σ < 0.0075 and by steepness
s t = k , l k 2 S k , l 1 / 2
in the range 0.05 < s t < 0.15 and by dispersion of Laplacian Λ = Δ η
σ L = Λ Λ ¯ 2 ¯ 1 / 2
in the range 6 < σ L < 12 . The calculations performed with such wide variations of integral parameters showed that the values A variate in the range 0.0010 < A < 0.0065 , whereas the dependence of w on w ζ is not linear as assumed in [11]. Hence, the problem is reduced to finding the dependence
A = w w ζ = f ( σ , σ L , s t )
In fact, several additional parameters were tested, whereas σ (Equation (27)) and σ L (Equation (29)) only turned out to be well connected with the ratio w / w ζ .
For generalization of the linear dependence (25) it is necessary to choose at least two parameters with a dimension of length taking into account internal characteristics of a wave field type (27)–(29). Finally, the best results were obtained in the following form
A = F ( μ ) ,
where μ is a parameter
μ = σ σ L ,
and function F is approximated by a formula
F = d 0 μ + d 1 μ + d 2 ,
where d 0 = 0.535 , d 1 = 0.0414 , d 2 = 0.00321 The function F μ is shown in Figure 2.
An approximation (33) is also correct for dimensional variables, because it is independent of the external scale L . The form of relation (31) and (32) and the constants in (33) were found on the basis of numerous numerical experiments with FWM (Equations (4)–(6)) and empirical selection of nondimensional variables, as well as functions and numerical parameters. That is why such simplified model is called the Heuristic Wave model (HWM).
The two-dimensional model, which is presumably able to replace the full 3D model, (4), (5), and (12), looks as follows:
η τ = η ξ φ ξ η ϑ φ ϑ + 1 + s w ,
φ τ = 1 2 φ ξ 2 + φ ϑ 2 1 + s w 2 η p
w ˜ = A 2 η ξ w ξ + η ϑ w ϑ + Δ η w s w ¯ ζ 1 + s .
Note that the right-hand side of Equation (36) contains full vertical velocity w = w ¯ + w ˜ as well as a linear component of vertical velocity w ¯ . An equation is represented in a form that is convenient for iterations. The iterations were performed until the next condition was reached.
max w ˜ i R i 1 < 10 7 ,
where R is the right-hand side of Equation (36), whereas i is the number of iterations. In the majority of cases, the number of iterations was equal to two and never exceeded four. Note that fulfilment of condition (37) can be regulated by an appropriate choice of parameters in the breaking algorithm (see below).
A comparison of function F calculated in the course of simulation by FWM with the results of calculations by Formula (33) is given in Figure 3. The agreement between these functions is quite satisfactory: a correlation coefficient between the model values F m and those calculated by Equations (31)–(33) is equal to 0.989, the root-mean-square difference between them equal to 0.014.
In fact, the research described in this paper was initiated due to the fact discovered in the first version of the model [11]. It was found that a nonlinear component of vertical velocity w ˜ (which is normally taken from solution of the Poisson equation) is a small value as compared with a full component w . However, the evolutionary Equations (34) and (35) contain a full component w = w ¯ + w ˜ only. It is illustrated in Figure 4 by comparison of the cumulative probability both for the full component w and the nonlinear components w ˜ (in Figure 4 they are marked as w t and w n , respectively). On the average, w ˜ is by one decimal order smaller than w . The main idea of the current approach is to replace the cumbersome algorithm for calculation of the vertical velocity equation φ ˜ with a simple 2D Equation (36). Such trick reduces the total 3D problem to the 2D one, which saves us from the endless problems associated with solving Poisson equations numerically, as well as from any other ways of reproduction of 3D structure of the velocity potential. It allows us to create a model running much faster than the 3D model.
The most convincing is the comparison of two independent methods of calculation of vertical velocity w , i.e., a method based on solution of Poisson Equation (12) (marked as w m ) and another one—based on the surface condition (36) and Equations (31)–(33) (marked as w c ). That comparison is given in Figure 5. The number of pairs of w m and w c used is 2 10 9 . The root-mean-square difference between w m and w c is equal to 0.0017.
The question is how well these simplified equations reproduce wave dynamics.
A simplified model is developed on the basis of the ‘exact’ model, both models having a similar structure and a set of parameters. The evolutionary equations of both models (4), (5), (34) and (35) are essentially the same. The difference between the models is determined by the simplified calculation of small nonlinear correction to the full vertical velocity introduced by Equations (31) and (33). A straightforward way of validation of the simplified model is performing the runs with the identical setting and initial conditions. A simple case of a quasi-stationary regime was considered.
The similarity of the models can be illustrated by difference D t (t is time)
D ( t )   =   σ f 2 + σ a 2 1 η i , j f η i , j a 2 ¯
and a correlation coefficient C t between the fields of elevation generated by a full model
C ( t )   =   σ f 2 σ a 2 1 / 2 η f η a ¯
where σ f and σ a are dispersions of the fields η f and η a ; the over-line denotes the averaging over the entire field. The right hand-sides of (38), (39) are calculated with the time interval equal to 100 Δ t . The evolution of D and C given in Figure 6 shows that the fields preserve their closeness well up to t = 100 (10,000 time steps). Such a result turned out to be quite unexpected because such a close agreement was observed for large wave fields composed of thousands of nonstationary and dispersing Fourier modes. The comparison of FWM with a simplified model was also made in [11], Figure 5. Both models used the spectral resolution of 129 × 129 Fourier modes and 512 × 256 knots. In the initial conditions JONSWAP spectrum with the peak wave number k p = 30 is assigned. It was obtained that the evolution of a correlation coefficient and a root-mean-square difference between the elevation fields was quite similar to those shown in Figure 6. Figure 6 proves that both models launched with the same boundary conditions produce a similar development for hundreds of wave peak periods. The advantage of a new scheme (Equations (31)–(33)) is its applicability to a wide range of external parameters.
The same comparison was made for FWM and HWM with a zero nonlinear correction to the vertical velocity w ˜ , i.e., using no Equation (36). The run was quite stable but an agreement between FWM and such fastest version of a simplified model was observed at a two- or three-times shorter period. Essentially, it is not clear what is more responsible for the nonlinearity: either the presence of nonlinear terms in the evolutionary Equations (34) and (35), or introducing a nonlinear correction to the vertical velocity through the exact Equation (12) or a simplified algorithm ((31)–(33) and (36)).

3. An Example of Simulation of the Wave Field Development on the Basis of HFM

The comparison of the long-term simulations performed with FWM and HWM was made in a previous publication [11]. That comparison proves that the accelerated and full models give the results quite similar to each other. The statistical characteristics of the solution were particularly well reproduced.
In the current work the calculations were performed with HWM model with the initiated algorithms for input and dissipation of energy. The algorithms for calculation of those terms are described in detail in [1,7,8,9]. The last minor modification of an algorithm for breaking is carried out in [10]. Those algorithms are not described in the current paper.
The model uses 257 × 257 Fourier modes or 1024 × 512 knots. The number of degrees of freedom (a minimum number of Fourier coefficients for elevation η and the surface velocity potential φ ) is 264,196. In the initial conditions the JONSWAP spectrum [12] with a peak wave number k p = 100 and the angle distribution proportional to sch ψ 256 was assigned. It corresponds to nearly unidirected waves with very small energy. Note that the details of the initial conditions assigned at sufficiently high wave numbers are not important, since the model quickly develops its own wave spectrum.
The model, as well as FWM, is supplied with a developed system of algorithms for monitoring the solution and recording many functions in a spectral and grid form. The most important were the spectra of elevation and different terms describing input and dissipation of energy. An exact expression for total kinetic energy of wave motion E k in the surface-following coordinates (1) can be represented in terms of surface variables φ and w = φ ζ :
E k = 0.5 1 + s φ w + 2 η ξ φ φ ξ + η ϑ φ φ ϑ + Δ η φ 2 ¯
where a single bar denotes the averaging over the ξ and ϑ coordinates. A relation (40) is easily derived with a use of Poisson Equation (7) and an expression for kinetic energy written in the coordinates (1). The total potential energy E p is simply calculated by the formula
E p = 0.5 η 2 ¯ ,
An equation of the integral energy E = E p + E k evolution can be represented in the following form:
d E d t = I ¯ + D b ¯ + D t ¯ + N ¯
where I ¯ is the integral input of energy from wind; D b ¯ is a rate of energy dissipation due to wave breaking; D t ¯ is a rate of energy dissipation due to filtration of the high-wave number modes (‘tail dissipation’); N ¯ is the integral effect of the nonlinear interactions described by the right-hand side of the equations when the surface pressure p is equal to zero. The differential forms for calculation of energy transformations can be, in principle, derived from Equations (4)–(6), but here a more convenient and simple method is applied. Different rates of the integral energy transformations can be calculated in the process of integration. For example, the value of I ¯ is calculated by the following relation:
I ¯ = 1 δ t E t + δ t E t ,
where E t + δ t is the integral energy of a wave field obtained after one time step with the right hand-side of Equation (35) containing only the surface pressure. This method is also applicable to the calculation of 1D spectra (i.e., the averaged over lateral wave number or along the direction in the polar coordinates) of different characteristics (see [1]).
An evolution of integral characteristics as a function of time t (expressed in the number of peak wave period t p = 2 π k p 1 / 2 in the initial conditions) is given in Figure 7. The fastest change of all the characteristics occurs in the first 400 t p until the moment when the breaking dissipation D b ¯ (curve 2) starts to develop. The input energy I ¯ (curve 1) remains nearly constant starting from t = 1000 . The total balance I ¯ + D b ¯ + D t ¯ (curve 4) gradually decreases but it does not turn into zero, and the total energy E ¯ ¯ (curve 4) continues to grow. Curve (5) describes an evolution of the spectrum-weighted wave number k s :
k s = k S d k d l S d k d l 1 / 2
It is a more convenient characteristic of the position of the spectrum because it changes in time smoothly. Over the time of integration, k s decreases from the value k s = 120 to k s = 20 .
The nonlinear transformation of spectrum described by the right hand-sides of Equations (34) and (35) always produces the output of energy outside the computational domain in Fourier space. This process can be treated as dissipation. In our opinion, to keep the whole process under control, it is reasonable to maintain the total energy in the adiabatic part of the problem, governed by the initial Equations (4)–(6) with no input energy and dissipation. It can be performed by correction of all Fourier coefficients for the potential f k , l and elevation h k , l :
f k , l = λ f k , l , h k , l = λ h k , l ,
where a coefficient λ is calculated by the formula
λ =   E f , h / E f , h 1 / 2 ,
where E is full energy calculated by a set of the coefficients f , h before a time step and by f , h —after a time step. The correction (45) and (46) is performed at each time step before calculating and adding the input and dissipation terms. A typical value of the coefficient λ is 1 ± 10 6 . It is important that the correction (45) and (46) eliminates not only the loss of energy due to energy outflow from the computational domain but also due to the errors of approximations. A typical value of integral dissipation of this type was by around 2 decimal orders smaller than other energy transformation effects. Finally, the integral effect of nonlinear interactions was equal to zero (a thin straight line in Figure 8).
Additional integral characteristics are shown in Figure 8 where the computed values as a function of fetch F . are compared with the same values calculated with a use of JONSWAP approximation [12]. Since the periodic variables are used, there is no other alternative but to assume that the fetch is a distance passed by a peak of the wave spectrum:
F ( t )   = 0 t c p d t
The dependence of F on t can be obtained with a use of JONSWAP approximations:
F = 3.34 · 10 3 t 3 / 2
The dependences (47) calculated by HWM and Eqauation (48) are shown in Figure 8 by thick and thin curves, respectively. As seen, these curves are quite close to each other. The dependence of peak frequency (defined as the frequency in a maximum of wave spectrum) is shown by curves 2. A dashed thick curve is calculated with the HWM model, and a thin curve corresponds to the dependence following from JONSWAP approximation:
ω p = 29.8 F 1 / 3
Opposite to Equation (44), the frequency ω p is a point value. Due to the irregularity of the 2D wave spectrum, the peak frequency changes in time not monotonically but in jumps, and sometimes ω p even moves back, i.e., to higher wave numbers. This is why the dependence ω F looks distorted. In Figure 8 this curve is smoothed. However, an agreement with the experimental dependence is quite obvious. The JONSWAP approximation of dependence of energy E p on fetch F is linear:
E p = 1.34 · 10 7 F
This dependence is shown in Figure 8 by a dotted line 3, whereas the dependence calculated by EWM is shown by thick dots. Those dependences agree with each other only at early stages of wave development. Evidently, an approximation (50) cannot be correct for a long fetch because it does not take into account the saturation of the spectrum. The more reasonable form is:
E p = f Ω p F ,
where f is a function of the inverse wave age Ω p = U / c p ( U is wind velocity).
In general, the results presented in Figure 7 and Figure 8 prove that an evolution of integral characteristics is reproduced satisfactorily.
The HWM also reproduces well the probability of the nondimensional surface elevation η n = η / H s . The height of large waves predicted by HWM is somewhat bigger than that predicted by FWM [9] due to a larger total steepness of wave field and because the integration time was much longer. It is interesting to note that normalization of surface elevation by the significant wave height is so efficient that the probability distribution for the nondimensional surface elevation P ( η n ) turns out to be nearly universal, i.e., it does not depend on wave energy and just slightly depends on total wave steepness. Curves in Figure 9 are calculated over 3.1 · 10 8 points; hence, it is not surprising that they are very smooth. In reality, the statistics of surface elevation shows very large scatter. The probability of the trough-to-crest height was calculated by a method of moving windows (see [13,14]). Each grey curve in Figure 10 is calculated over individual wave field containing 524.288 points. The total number of curves is 594. As seen, the probability for different wave fields variates in a broad range: many of them do not contain ‘freak’ waves with the height of η t c > 2 at all, whereas in many cases their height can exceed such a high value as η t c = 2.5 and even reach a value η t c = 3 . This is why prediction of extreme waves for practical purposes is senseless with no indication of confidence intervals or distribution of probability for each value of η t c .
Figure 11 demonstrates statistical characteristics of positive and negative elevations as a function of the trough-to-crest height. As seen, the crest height Z c grows with the increase in total height Z t c faster than the trough depth, which can be explained by nonlinearity of the process. Figure 11 shows that on the average the trough-to-crest height Z t c / H s = 2 corresponds to the crest height Z c / H s = 1.3 and the trough depth Z t = 0.7 . The dependences between those characteristics can be approximated by the formulas:
Z ˜ c = 0.006 + 0.4162 Z ˜ t c + 0.111 Z ˜ t c 2 Z ˜ t = 0.006 0.5838 Z ˜ t c + 0.111 Z ˜ t c 2
where tilde denotes normalization by the significant wave height H s . Note that the first and third coefficients in approximation (52) coincide, which proves high quality of the ensemble. A dashed curve corresponds to the same data obtained with FWM [10]. The difference between the new and old data is growing for large values of Z ˜ c . The old data contain 41.677 points, whereas a new array contains 1.264.441 points.
In a linear wave field (i.e., superposition of linear modes with random phases), the probabilities of crest height and trough depths are identical, but the probability of the trough-to-crest height is close to that for the nonlinear waves in a wave field with the same energy. This fact reliably proved with extensive calculations using the accurate two-dimensional and three-dimensional models [13,14,15], has not been explained yet. This fact seems to be directly related to the physics of extreme waves. The nonlinearity manifests itself in vertical asymmetry of waves, but not in the probability of full height.
The next three figures show the 2D spectra obtained by averaging over six consecutive periods of length t = 1000 (about 1600 initial peak period). Each spectrum is calculated by averaging 100 spectra transformed into the polar coordinates θ , k . The surface elevation spectra S θ , k are given in Figure 12. The colors correspond to different values of log 10 S / S m where S m is a maximum of the last of the six spectra. A minimum value (black color) corresponds to the values less than 10 4 S m . It is seen that the spectrum grows and approaches the point k = 0 . Despite the severe averaging, the spectra look highly irregular. They contain multiple peaks and holes (well seen in panels 5 and 6). The discussion of this effect is given in [7,13].
The spectra of the energy input shown in Figure 13 calculated similar to wave spectrum are much wider than the wave spectrum, because the input of energy decreases with growth of wave number slower than the wave spectrum. It happens because β -function (see [16]) depends on wave number quadratically. For large fetches the areas of negative wind input in low wave numbers (i.e., a flux of energy from fast wave modes to wind) usually arise, but it cannot be shown in the given figures.
The spectra of total dissipation rate equal to the sum of breaking and tail dissipation, are shown in Figure 14. A maximum of dissipation always falls on a maximum of wave spectrum. However, dissipation, as well as the input energy, occurs everywhere including the spectral tail. Here, as in all our previous calculations, we did not reproduce an adiabatic regime at the spectral tail, which excludes any possibility to draw analogy to Kolmogorov’s spectrum. The attenuation of energy in the spectral tail depends on a shape of entire spectrum, as well as on wind velocity, wave age and other factors.
The process of downshifting is well seen in Figure 15 where the evolution of 1D wave spectra (black curves) and the spectra of the nonlinear transformation rate (red curves) in the vicinity of a peak averaged over three consequent periods, each of t = 2000 length, are given. All the spectra are calculated over 198 records of the two-dimensional spectra transformed into the polar coordinates, summated over angles and represented as a function of radius r = k . The locations of a wave spectrum maximum are indicated by thin vertical lines. Note that the analysis of the spectra for a nonlinear transformation rate is not simple, because the dispersion of reversible exchanges of energy between the modes is much higher than the rate of residual irreversible interaction. The accuracy of calculations grows with increase in spectral resolution. The recently launched calculations with HWF using spectral resolution of 1027 × 1027 modes gave the averaged spectral distribution much smoother than that in Figure 15. In general, the spectra in Figure 15 agree with the known regularities: the energy from a right slope of wave spectrum is transferred to the left slope, i.e., to longer waves. It qualitatively confirms Hasselmann’s theory [17] considering the four-wave interactions. Opposite to that theory, actual interactions occur between all the spectral modes within the scope of the Fourier transform method. Unfortunately, the quantitative comparison with the theory is impossible since neither of the codes for calculation of Hasselmann integral can operate with such a high number of modes as used in phase-resolving models.
Most of the energy coming from wind to waves dissipates because of wave breaking. A currently used parameterization method for wave breaking is based on a diffusion algorithm with a variable coefficient of diffusion applied to both the elevation and surface velocity potential. The breaking occurs when Laplacian of elevation Δ = η ξ ξ + η ϑ ϑ becomes smaller than the negative value Δ c (in the current model Δ c = −80). A diffusion coefficient B is defined as follows:
B = C b Δ x 2 Δ 2 Δ < Δ c 0 Δ Δ c ,
where Δ x is a horizontal step, and C s = 0.02 is a coefficient. The parameters Δ c and C b are selected with consideration of a rate of growth of total energy. It is remarkable that the integral rate of breaking dissipation is not too sensitive to the choice of parameters because the effects of frequency of breaking (regulated by Δ c ) compensate for the intensity of breaking (depending on C b ). The rate of breaking seems to depend mostly on a rate of the energy input: most of the energy should be removed to maintain the numerical stability in the model, i.e., to prevent appearance of too sharp waves. Obviously, in a real wave field similar effects occur: dissipation absorbs extra energy obtained from the wind. Though a suggested algorithm seems very simple, the scheme performs all that is required.
Note that some part of the released horizontal momentum and energy should be transferred to horizontal motion (current). This effect cannot be taken into account in a potential model. The energy released after breaking is partially dissipated and partially distributed in space and between the wave modes, thus additionally contributing to the nonlinear interactions. The effects of breaking are presented in Figure 16 where the number of breaking cases (% of the total number of knots) and the integral rate of energy loss at each time step are presented. As seen, both characteristics grow in time with the tendency for stabilization by the end of a run.
At an early stage of the work on parameterization of breaking the attempts were made to construct a scheme based on the geometrical and mechanical consideration with direct calculations of the stability criterion, a detaching volume of water and the possible consequences. Finally, we came to the conclusion that such a scheme would never work properly or be sufficiently robust to be used at millions of grid point and time steps. It is more reliable to adhere to the standard methods adopted in the geophysical fluid mechanics.
The width of spectrum can be characterized by a function Ψ ω / ω p (see [18,19]).
Ψ = S ω , ψ ψ d ω d ψ S ω , ψ d ω d ψ
where integrals have taken over the domain 0 < ω < ω c , π / 2 < ψ < π / 2 . The wave spectra as the functions of frequency ω normalized by peak frequency ω p for the same three periods shown in Figure 15 are given. As seen, the Ψ curves corresponding to different wave ages are close to each other. All of them have a sharp maximum at the frequencies below the spectral peak, as well as a well-pronounced minimum in the spectral peak and a relatively slow growth above the spectral peak. The decrease in Ψ at high frequencies is probably caused by the high-frequency dumping. The angle distribution was investigated in [18,19]. The approximations of Ψ ω / ω p from the different sources collected in [20] show considerable scatter, but the general features are quite similar to those shown in Figure 17.
An example of a simulated elevation field η normalized by the significant wave height H s is given in Figure 18. In the initial conditions, small unidirected waves at inverse wave age U / c p were assigned. In Figure 18 we can see a developed three-dimensional field with the energy by two decimal orders higher than that in the initial conditions, and the broad angle distribution. Black dots mark the points where breaking occurs at a given time step.

4. Conclusions

The three-dimensional modeling of surface waves based on full nonlinear equations is a powerful tool for the investigation of wave processes, development of the parameterization schemes for the phase-resolving and spectral models and direct simulation of wave regimes in small basins. This type of modeling is rapidly developing. However, all approaches have a common limitation, i.e., low computational efficiency. Working with such models, even with modest resolution, turns into an endless waiting of the results. This property of the models slows down their improvement, in particular, development of the parameterization schemes for physical processes, since such work requires multiple repetitions of the runs. Such schemes significantly depend on resolution of a model, which limits the possibility of the research with low-resolution models.
The current paper continues to develop a new approach to the phase-resolving modeling of two-dimensional periodic wave fields at infinite depth. The basic concept of the scheme follows from presentation of the velocity potential as a sum of the linear and nonlinear components suggested in [7]. The solution for a linear component is known; hence, a nonlinear component should be calculated through Poisson equation with a zero-boundary condition on the surface. It was observed in [10] that such approach offers a new way to simplify the calculation by considering of 2D Poisson equation on the surface. The equation that can be treated as additional exact surface condition contains both the first φ ζ and second φ ζ ζ vertical derivatives of the potential. Thus, the system of equations remains unclosed. It was empirically discovered that these variables are closely connected to each other. The linear dependence between φ ζ and φ ζ ζ was checked in [10]. It was shown that the use of such hypotheses leads to formulation of a closed system of equations.
In this paper, the above approach was improved by introducing a more flexible dependence of a ratio φ ζ / φ ζ ζ on the nondimensional integral parameters of the wave field. The comparison of full vertical velocity w = φ ζ calculated using a new method, with that calculated by FWM, shows good agreement (Figure 5). The point-to-point comparison of runs with HWM and FWM starting from the identical initial conditions proves that the solutions agree over thousands of time steps (Figure 6). A long high-resolution run with HWM revealed that a simplified model produces the statistical results similar to those obtained with FWM [9]. The main advantage of the new approach is high computational efficiency. The ratio of time spent for calculation with identical FWM and HWM variated in the range of 50–100. This ratio depends on complexity of parameterization of physics and volume/frequency of the processed and recorded information. Some acceleration can be reached by replacing a Runge–Kutta time integration scheme with a high accuracy one-step scheme.
The model suggested is not mathematically precise since it is based on the empirical (in computational sense) connection between the local variables, i.e., vertical velocity and its vertical derivative. The fact that those variables are connected by an equation containing horizontal derivatives indicates that the local connection cannot be exact. To discuss this problem, it is necessary to answer a general question: what does the concept of precision mean?
It should be noted that surface waves are often investigated on the basis of different simplified equations such as a linearized Navier–Stokes equation, Korteweg de Vrie equation, a linear and nonlinear Schrödinger (NLS) equation and others. No one actually knows which real properties of sea waves are omitted in those models. In some simplified approaches (e.g., in the models based on a nonlinear Schrödinger equation (see [21]), the instability of high waves is missing; hence, the amplitudes of simulated waves can be unrealistically large. In reality, an excessive growth of waves is prevented due to the energy focusing and wave-breaking instability. The models based on the full potential equations predict the probability distribution for large waves quite satisfactorily. The main obstacle to a widespread use of the models based on complete equations is their low computational efficiency. An approach suggested here eliminates this obstacle by excluding reproduction of a three-dimensional flow structure. However, nothing prevents us from using a combined model with the on-demand calculations of a three-dimensional field of the velocity potential. A simplified model gives almost the same statistical results as a full model. Essentially, the difference between the statistical results given by 3D and 2D models are similar to those obtained with a 3D model initiated at different sets of random phases. It is easy to see that the equations suggested above are completely similar to the full equations, with one exception, i.e., a small non-linear correction to the total vertical velocity is calculated not from a 3D Poisson Equation (12) but from a simple 2D Equation (17). It can be added that a simplified model has a much simpler structure than FWM and is easily programmed. It is known that in most cases the data on 3D velocity potential are not required, but, if necessary, they can be easily obtained periodically in a process of model run or by the post-processing using the restart records and the corresponding subroutines adopted from 3D FWM.
It is obvious that an approach developed here is not fully precise. For example, it cannot be applied for individual cases with a small number of modes, for example, for simulation of a steep Stokes wave, as it was demonstrated for FWM in [7,9]. The 2D model is intended for simulation of the statistically homogeneous multi-mode wave field.
To continue the discussion on the accuracy problem, it should be noted that the full 3D adiabatic model is also inaccurate, because it is suitable for relatively short time intervals only. Real waves receive energy from wind and dissipate. Currently, such processes are poorly studied because they are evidently more complicated than the wave movement itself. Therefore, neither of the models can be considered as fully adequate. Especially for the followers of high accuracy it should be added that the full potential equations are far from being accurate even for the waves with no input energy and breaking. Real waves develop in a flow with the velocity shift, and they exchange energy with vortex movements [22]. Therefore, they are better described by the Navier–Stokes equations. As a result, we obtain a variety of problems associated with the turbulence. It may turn out that the simplifications used in the construction of this model may be quite insignificant, as compared with other disadvantages.
Certainly, there are many ways to improve the model, first of all, for example, testing of alternative schemes for calculation of a nonlinear velocity component (Equations (31)–(33)). The attempts were made to construct the closure scheme for a surface Poisson equation in Fourier space. This way turned out to be unrealistic, since the local structure of the velocity potential depends not only on Fourier variables but also on a local full elevation.
The most evident advantage of a new model is the absence of calculation of the 3D structure of the velocity potential, which increases the speed of calculations by around two orders. This easily programmable model can be used to improve the physical components of the model and the long-term simulations of a multimode wave field. The model can be also used for acceleration of the calculations as a block of an exact 3D model (such as FWM). This coupling provides the possibility of runtime tuning of the parameters of the closure scheme for a surface Poisson equation.
The scope of practical applications of such model is quite wide. First of all, it should be used for building up of parametrization of physical processes in spectral and phase-resolving models. HWM, due to its high performance, can be used as a component of wave prediction models (such as WAM or WAVEWATCH) for interpretation of spectral information in terms of the phase-resolved wave field. A possible consecution of action is as follows: the spectra calculated in the spectral model in the selected areas are used for generation of the initial conditions in the wave number space; then the calculations up to the stationary statistical regime with HWM are performed; finally, the statistical characteristics of wave field including the probability distribution of wave height are calculated. The input and dissipation terms can be included in such simulation. Note that such approach allows calculating nonlinear interactions between all the modes based on full nonlinear equations. Such scheme can significantly extend the output of spectral modeling.
It is comparatively easy to transform the scheme into a finite-difference form and use the model for investigation of a wave regime in a specific geographical environment with variable boundary conditions. Note that a numerical scheme for this model is much closer to the local scheme than FWM, which makes convenient to carry out the parallel computations. An investigation of applicability of a developed approach for the waves above finite depth is underway.

Funding

This research was supported by Russin Science Foundation RSF (project No. 22-21-00139).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The author would like to thank O. Chalikova for her assistance in preparation of the manuscript.

Conflicts of Interest

The author declares no conflict of interest.

References

  1. Chalikov, D. Numerical modeling of surface wave development under the action of wind. Ocean Sci. 2018, 14, 453–470. [Google Scholar] [CrossRef] [Green Version]
  2. Grilli, S.; Guyenne, P.; Dias, F. A fully non-linear model for three-dimensional overturning waves over an arbitrary bottom. Int. J. Numer. Meth. Fluids 2001, 35, 829–867. [Google Scholar] [CrossRef]
  3. Xu, L.; Guyenne, P. Numerical simulation of three-dimensional nonlinear water waves. J. Comput. Phys. 2009, 228, 8446–8466. [Google Scholar] [CrossRef]
  4. Grue, J.; Fructus, D. Model for Fully Nonlinear Ocean Wave Simulations Derived Using Fourier Inversion of Integral Equations in 3D. In Advances in Numerical Simulation of Nonlinear Water Waves; University of London: London, UK, 2010; pp. 1–42. [Google Scholar]
  5. Engsig-Karup, A.; Madsen, M.; Glimberg, S. A massively parallel GPU-accelerated mode for analysis of fully nonlinear free surface waves. Int. J. Numer. Methods Fluids 2012, 70, 20–36. [Google Scholar] [CrossRef] [Green Version]
  6. Ducrozet, G.; Bonnefoy, F.; Le Touzé, D.; Ferrant, P. HOS-ocean: Open-source solver for nonlinear waves in open ocean based on High-Order Spectral method. Comp. Phys. Commun. 2016, 203, 245–254. [Google Scholar] [CrossRef] [Green Version]
  7. Chalikov, D.; Babanin, A.V.; Sanina, E. Numerical Modeling of Three-Dimensional Fully Nonlinear Potential Periodic Waves. Ocean. Dyn. 2014, 64, 1469–1486. [Google Scholar] [CrossRef]
  8. Chalikov, D. Numerical Modeling of Sea Waves; Springer: Cham, Switzerland, 2016; p. 330. [Google Scholar] [CrossRef]
  9. Chalikov, D. High-Resolution Numerical Simulation of Surface Wave Development under the Action of Wind. In Geophysics and Ocean Waves Studies; IntechOpen: London, UK, 2020. [Google Scholar] [CrossRef]
  10. Thomas, L.H. Elliptic Problems in Linear Differential Equations over a Network; Columbia University: New York, NY, USA, 1949. [Google Scholar]
  11. Chalikov, D. Accelerated reproduction of 2-D periodic waves. Ocean. Dyn. 2021, 71, 309–322. [Google Scholar] [CrossRef]
  12. Hasselmann, K.; Barnett, T.P.; Bouws, E.; Carlson, H.; Cartwright, D.E.; Enke, K.; Ewing, J.A.; Gienapp, A.; Hasselmann, D.E.; Kruseman, P.; et al. Measurements of wind-wave growth and swell decay during the Joint Sea Wave Project (JONSWAP). Ergaenzungsheft zur Deutschen Hydrographischen Zeitschrift Reihe A 1973, A8, 1–95. [Google Scholar]
  13. Sanina, E. Statistics of Wave Kinematics in Random Directional Wave Fields. Bull. Aust. Math. Soc. 2015, 93, 169. Available online: http://hdl.handle.net/1959.3/394825 (accessed on 12 February 2022). [CrossRef] [Green Version]
  14. Chalikov, D. Freak waves: Their occurrence and probability. Phys. Fluid 2009, 21, 076602. [Google Scholar] [CrossRef]
  15. Chalikov, D.; Babanin, A.V. Comparison of linear and nonlinear extreme ave statistics. Acta Oceanol. Sin. 2016, 35, 99–105. [Google Scholar] [CrossRef]
  16. Chalikov, D.; Rainchik, S. Coupled numerical modelling of wind and waves and the theory of the wave boundary layer. Bound.-Layer Meteorol. 2010, 138, 1–41. [Google Scholar] [CrossRef]
  17. Hasselmann, K. On the non-linear energy transfer in a gravity wave spectrum, part 1. J. Fluid Mech. 1962, 12, 481–500. [Google Scholar] [CrossRef]
  18. Donelan, M.A.; Hamilton, J.; Hui, W.H. Directional spectra of wind-generated waves. Philos. Trans. R. Soc. 1985, 315, 509–562. [Google Scholar] [CrossRef]
  19. Babanin, A.V.; Soloviev, Y.P. Field investigation of transformation of the wind wave frequency spectrum with fetch and the stage of development. J. Phys. Oceanogr. 1998, 28, 563–576. [Google Scholar] [CrossRef]
  20. Liu, Q.; Rogers, W.E.; Babanin, A.V.; Young, I.R.; Romero, L.; Zieger, S.; Qiao, F.; Changlong, G. Observation-Based Source Terms in the Third-Generation Wave Model WAVEWATCH III: Updates and Verification. J. Phys. Oceanogr. 2019, 49, 489–517. [Google Scholar] [CrossRef]
  21. Slunyaev, A.; Kharif, C.; Pelinovsky, E.; Talipova, T. Nonlinear wave focusing on water of finite depth. Physics D 2002, 173, 77–96. [Google Scholar] [CrossRef]
  22. Babanin, A.; Chalikov, D. Numerical Investigation of Turbulence Generation in Non-Breaking Potential waves. J. Geophys. Res. 2012, 117, C00J17. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Dependence of a nonlinear component of vertical velocity w ˜ (vertical axis) on a vertical derivative φ ˜ ζ ζ = w ˜ ζ . The thick line is the averaged over all points values; thin lines show dispersion. Both are calculated by averaging by bins of size 0.01. The dashed line shows the probability distribution for w ζ normalized by a maximum of probability. Value w ˜ ζ = 0.02 corresponds to the zero probability; w ˜ ζ = 0.02 corresponds to the probability equal to 1.
Figure 1. Dependence of a nonlinear component of vertical velocity w ˜ (vertical axis) on a vertical derivative φ ˜ ζ ζ = w ˜ ζ . The thick line is the averaged over all points values; thin lines show dispersion. Both are calculated by averaging by bins of size 0.01. The dashed line shows the probability distribution for w ζ normalized by a maximum of probability. Value w ˜ ζ = 0.02 corresponds to the zero probability; w ˜ ζ = 0.02 corresponds to the probability equal to 1.
Jmse 10 00410 g001
Figure 2. Dependence of a nondimensional variable F = A / σ on a nondimensional parameter μ = σ σ L (grey points). Solid line corresponds to the approximation (32).
Figure 2. Dependence of a nondimensional variable F = A / σ on a nondimensional parameter μ = σ σ L (grey points). Solid line corresponds to the approximation (32).
Jmse 10 00410 g002
Figure 3. Comparison of function F m calculated in the course of simulation by FWM with the results of calculation of F c by Equations (31)–(33) (grey dots). The thick line is calculated by averaging by bins of size 0.01.
Figure 3. Comparison of function F m calculated in the course of simulation by FWM with the results of calculation of F c by Equations (31)–(33) (grey dots). The thick line is calculated by averaging by bins of size 0.01.
Jmse 10 00410 g003
Figure 4. Cumulative probability of a full vertical velocity component w (marked by w t ) and a nonlinear component of vertical velocity w ˜ ( w n ).
Figure 4. Cumulative probability of a full vertical velocity component w (marked by w t ) and a nonlinear component of vertical velocity w ˜ ( w n ).
Jmse 10 00410 g004
Figure 5. Comparison of full vertical velocity w m calculated by FWM with the same velocity w c calculated by Equations (31)–(33). Thick line corresponds to the averaged dependence, whereas thin lines show the dispersion. Both are calculated by averaging by bins of size 0.001. Dashed curve shows the probability distribution for w normalized by a maximum value of the probability.
Figure 5. Comparison of full vertical velocity w m calculated by FWM with the same velocity w c calculated by Equations (31)–(33). Thick line corresponds to the averaged dependence, whereas thin lines show the dispersion. Both are calculated by averaging by bins of size 0.001. Dashed curve shows the probability distribution for w normalized by a maximum value of the probability.
Jmse 10 00410 g005
Figure 6. Evolution of a correlation coefficient C t and the root-mean-square difference D t between the elevation fields η calculated by HWM and FWM with the same initial conditions.
Figure 6. Evolution of a correlation coefficient C t and the root-mean-square difference D t between the elevation fields η calculated by HWM and FWM with the same initial conditions.
Jmse 10 00410 g006
Figure 7. An evolution of integral characteristics): 1—a rate of the input energy 10 7 I ¯ ); 2—a rate of energy loss due to the breaking 10 7 D b ¯ ; 3—a rate of energy loss due to the tail dissipation 10 7 D b ¯ ; 4—the balance of energy 10 7 I ¯ + D b ¯ + D t ¯ ; 5 is the weighted peak wave number (multiplied by 10 3 ); Equation (44); 6—the potential energy (multiplied by 10 7 ). A thin straight line demonstrates a zero integral effect of the nonlinear interactions.
Figure 7. An evolution of integral characteristics): 1—a rate of the input energy 10 7 I ¯ ); 2—a rate of energy loss due to the breaking 10 7 D b ¯ ; 3—a rate of energy loss due to the tail dissipation 10 7 D b ¯ ; 4—the balance of energy 10 7 I ¯ + D b ¯ + D t ¯ ; 5 is the weighted peak wave number (multiplied by 10 3 ); Equation (44); 6—the potential energy (multiplied by 10 7 ). A thin straight line demonstrates a zero integral effect of the nonlinear interactions.
Jmse 10 00410 g007
Figure 8. An evolution of integral characteristics calculated with HWM (thick curves) and JONSWAP approximation [12]: 1—dependence of fetch 10 4 F (horizontal axis, Equation (47)) expressed in terms of the initial peak wave length, on time 10 4 t (vertical axis, Equation (48)) expressed in the periods of the initial peak wave period; 2—dependence of the peak frequency 10 2 ω p on fetch 10 4 F ’ 3—dependence of the potential wave energy 10 5 E p on fetch 10 4 F .
Figure 8. An evolution of integral characteristics calculated with HWM (thick curves) and JONSWAP approximation [12]: 1—dependence of fetch 10 4 F (horizontal axis, Equation (47)) expressed in terms of the initial peak wave length, on time 10 4 t (vertical axis, Equation (48)) expressed in the periods of the initial peak wave period; 2—dependence of the peak frequency 10 2 ω p on fetch 10 4 F ’ 3—dependence of the potential wave energy 10 5 E p on fetch 10 4 F .
Jmse 10 00410 g008
Figure 9. Distribution of the probability of the nondimensional surface elevation calculated by HWM (thick curve) and FWM (thin curve).
Figure 9. Distribution of the probability of the nondimensional surface elevation calculated by HWM (thick curve) and FWM (thin curve).
Jmse 10 00410 g009
Figure 10. The cumulative probability of the nondimensional wave trough-to-crest height Z t c / H s . Grey curves are calculated for each individual 2D wave field with the dimensions of 1024 × 512 knots. Thick curve is the averaged probability over all 594 fields.
Figure 10. The cumulative probability of the nondimensional wave trough-to-crest height Z t c / H s . Grey curves are calculated for each individual 2D wave field with the dimensions of 1024 × 512 knots. Thick curve is the averaged probability over all 594 fields.
Jmse 10 00410 g010
Figure 11. Dependence of the crest height Z c / H s (top section) and the trough depth Z t / H s (bottom section) on the total wave height Z t c / H s . Solid curves show the quadratic approximation of the point data, whereas a dashed curve represents the same approximations obtained in [9].
Figure 11. Dependence of the crest height Z c / H s (top section) and the trough depth Z t / H s (bottom section) on the total wave height Z t c / H s . Solid curves show the quadratic approximation of the point data, whereas a dashed curve represents the same approximations obtained in [9].
Jmse 10 00410 g011
Figure 12. The 2D wave spectra averaged over six consecutive periods, each with length t = 2000 . All the spectra are normalized by a maximum of the last spectrum. Black color corresponds to zero.
Figure 12. The 2D wave spectra averaged over six consecutive periods, each with length t = 2000 . All the spectra are normalized by a maximum of the last spectrum. Black color corresponds to zero.
Jmse 10 00410 g012
Figure 13. The 2D spectra of the energy input averaged over six consecutive periods, each with length t = 2000 . All the spectra are normalized by a maximum of the last spectrum. Black color corresponds to zero.
Figure 13. The 2D spectra of the energy input averaged over six consecutive periods, each with length t = 2000 . All the spectra are normalized by a maximum of the last spectrum. Black color corresponds to zero.
Jmse 10 00410 g013
Figure 14. The 2D spectra of total dissipation averaged over six consecutive periods, each with length t = 2000 . All the spectra are normalized by a maximum of the last spectrum. Black color corresponds to zero.
Figure 14. The 2D spectra of total dissipation averaged over six consecutive periods, each with length t = 2000 . All the spectra are normalized by a maximum of the last spectrum. Black color corresponds to zero.
Jmse 10 00410 g014
Figure 15. The 1D wave spectra multiplied by 10 5 (black curves) and the spectra of a nonlinear transformation rate multiplied by 10 8 (red curves) in the vicinity of the inverse wave age: 1— U / c p = 1.50 ; 2— U / c p = 1.38 ; 3— U / c p = 1.20 .
Figure 15. The 1D wave spectra multiplied by 10 5 (black curves) and the spectra of a nonlinear transformation rate multiplied by 10 8 (red curves) in the vicinity of the inverse wave age: 1— U / c p = 1.50 ; 2— U / c p = 1.38 ; 3— U / c p = 1.20 .
Jmse 10 00410 g015
Figure 16. Thin curve is an evolution of a number of breaking cases at each time step (% of the total number of knots 524.288 in the elevation field); a dotted curve is a rate of the integral energy lost (with a reverse sign) due to the breaking (multiplied by 10 6 ).
Figure 16. Thin curve is an evolution of a number of breaking cases at each time step (% of the total number of knots 524.288 in the elevation field); a dotted curve is a rate of the integral energy lost (with a reverse sign) due to the breaking (multiplied by 10 6 ).
Jmse 10 00410 g016
Figure 17. The wave spectra as functions of the nondimensional frequency ω / ω p , averaged over three consecutive periods, each of t = 2000 length (black curves). The corresponding functions ϒ ω / ω p (Equation (53); ω p is a frequency in the spectral peak) (red curves).
Figure 17. The wave spectra as functions of the nondimensional frequency ω / ω p , averaged over three consecutive periods, each of t = 2000 length (black curves). The corresponding functions ϒ ω / ω p (Equation (53); ω p is a frequency in the spectral peak) (red curves).
Jmse 10 00410 g017
Figure 18. An example of an elevation field η normalized by the significant wave height H s . Black dots mark the points of wave breaking. A maximum of the normalized trough-to crest height in this field is equal to 2.1. Horizontal axis corresponds to ξ and vertical axis—to ϑ .
Figure 18. An example of an elevation field η normalized by the significant wave height H s . Black dots mark the points of wave breaking. A maximum of the normalized trough-to crest height in this field is equal to 2.1. Horizontal axis corresponds to ξ and vertical axis—to ϑ .
Jmse 10 00410 g018
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Chalikov, D. A 2D Model for 3D Periodic Deep-Water Waves. J. Mar. Sci. Eng. 2022, 10, 410. https://doi.org/10.3390/jmse10030410

AMA Style

Chalikov D. A 2D Model for 3D Periodic Deep-Water Waves. Journal of Marine Science and Engineering. 2022; 10(3):410. https://doi.org/10.3390/jmse10030410

Chicago/Turabian Style

Chalikov, Dmitry. 2022. "A 2D Model for 3D Periodic Deep-Water Waves" Journal of Marine Science and Engineering 10, no. 3: 410. https://doi.org/10.3390/jmse10030410

APA Style

Chalikov, D. (2022). A 2D Model for 3D Periodic Deep-Water Waves. Journal of Marine Science and Engineering, 10(3), 410. https://doi.org/10.3390/jmse10030410

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop