**5. Evolution of wave field**

The integration was done for 1,200,000 steps with the time step Δ*t* ¼ 0*:*005 up to the nondimensional time *T* ¼ 6000, which corresponded to 9550 initial wave peak periods. The total energy of wave motion *Е* ¼ *Ep* þ *Ek* (*Ep* is the potential energy, while *Ek* is the kinetic energy) is calculated with the following formulas:

$$E\_p = 0.25\overline{\eta^2}, \quad E\_k = 0.5\overline{\overline{\left(\rho\_x^2 + \rho\_y^2 + \rho\_z^2\right)}}, \quad E = E\_p + E\_k,\tag{27}$$

where a single bar denotes the averaging over the *ξ* and *ϑ* coordinates, while a double bar denotes the averaging over the entire volume. The derivatives in Ref.

*High-Resolution Numerical Simulation of Surface Wave Development under the Action of Wind DOI: http://dx.doi.org/10.5772/intechopen.92262*

(27) are calculated according to the transformation (1). An equation of the integral energy *E* evolution can be represented in the following form:

$$\frac{dE}{dt} = \overline{I} + \overline{D\_b} + \overline{D\_t} + \overline{N},\tag{28}$$

where *I* is the integral input of energy from wind (Eqs. (14)–(20)); *Db* is a rate of the energy dissipation due to wave breaking (Eqs. (23)–(26)); *Dt* is a rate of the energy dissipation due to filtration of high-wave number modes ('tail dissipation,' Eqs. (21) and (22)); *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 the energy transformations can be, in principle, derived from Eqs. (4)–(6), but here a more convenient and simple method was applied. Different rates of the integral energy transformations can be calculated with the help of fictitious time steps (i.e., apart from the basic calculations). For example, the value of *I* is calculated by the following relation:

$$\overline{\overline{I}} = \frac{1}{\Delta t} \left( \overline{\overline{E^{t+\Delta t}}} - \overline{\overline{E^t}} \right), \tag{29}$$

where *Et*þΔ*<sup>t</sup>* is the integral energy of a wave field obtained after one time step with the right side of Eq. (6) containing only the surface pressure calculated with Eqs. (14)–(18).

The evolution of the characteristics calculated by formula (29) is shown in **Figure 2**.

The sharp variation of all the characteristics at *t*<500 is explained by adjustment of the linear initial fields to the nonlinearity. The integral effect of the nonlinear interaction *I* (straight line 1) was very close to zero. The tail dissipation *Dt* (curve 2) is smaller than the breaking dissipation *Dt*(curve 3). The value of *Db* has significant fluctuations due to introduction of the criteria (25) and (26). The dissipation *Db* þ *Dt* absorbs nearly all of the incoming energy, and just a small part of it is going for growth of waves. The balance of energy *B* ¼ *I* þ *Dt* þ *Db* (curve 5) fluctuates and approaches zero when energy *E* (curve 6 in **Figure 2**) approaches saturation.

The time evolution of the integral spectral characteristics is presented in **Figure 3**.

Curve 1 corresponds to the weighted frequency *ω<sup>w</sup>*

$$
\rho\_w = \left(\frac{\int \rho \text{Sdd} \, dl}{\int \text{Sdd} \, dl}\right)^{1/2},\tag{30}
$$

where integrals are taken over the entire Fourier domain. The value *ω<sup>w</sup>* is not sensitive to the details of spectrum; hence, it well characterizes the position of spectrum and spectral peak shifting. Curve 2 describes the evolution of the spectral maximum. The step shape of curve corresponds to the fundamental property of downshifting. Opposite to common views, the development of spectrum occurs not monotonically, but by appearance of a new maximum at a lower wave number as well as by attenuation of the previous maximum. It is interesting to note that the same phenomenon is also observed in the spectral model [36].

The value of fetch in the periodic problem can be calculated by integration of the peak phase velocity *cp* <sup>¼</sup> j j *<sup>k</sup>* �1*=*<sup>2</sup> over time.

### **Figure 2.**

*The evolution of integral characteristics. The rate of evolution of the integral energy multiplied by* 10<sup>8</sup> *due to: (1) nonlinear interaction I (Eq. (29)); (2) tail dissipation Dt(Eqs. (21) and (22)); (3) breaking dissipation Db (Eqs. (23)–(26)); (4) input of energy from wind I (Eqs. (14)–(20)); and (5) balance of energy I* þ *Dt* þ *Db. Curve 6 shows the evolution of wave energy* 105*E. Vertical gray bars show instantaneous values; thick curve shows the smoothed behavior.*

$$F = \int\_{t\_0}^{t} c\_p dt. \tag{31}$$

The numerical experiment reproduces the case when development of wave field occurs under the action of a permanent and uniform wind. This case corresponds to the JONSWAP experiment [29]. It is suggested that the frequency of spectral peak changes as *F*�1*=*<sup>3</sup> , while the full energy grows linearly with *F:* Neither of the dependences can be exact since they do not take into account approaching a stationary regime. Besides, the dependence of frequency on fetch is singular at *F* ¼ 0*:* A more accurate is the approximation:

$$
\rho\_p = \frac{75.6}{5.63 + F^{1/3}}.\tag{32}
$$

Obviously, the dependence *<sup>ω</sup><sup>p</sup>* � *<sup>F</sup>*�1*=*<sup>3</sup> is valid in a narrower interval of *<sup>F</sup>*. As seen, contrary to *ωw*, the peak frequency changes not monotonically, but by appearance of a new maximum at a lower wave number as well as by attenuation of the previous maximum. It is interesting to note that the same phenomenon is also observed in the spectral model (16). The dependence of the total energy *E* on fetch *F* does not look like a linear one, but it is worth to note that the JONSWAP dependence is evidently inapplicable to a very small and large fetch.

*High-Resolution Numerical Simulation of Surface Wave Development under the Action of Wind DOI: http://dx.doi.org/10.5772/intechopen.92262*

**Figure 3.**

*Dependence of integral characteristics on fetch (Eq. (32)): (1) weighted mean frequency ω<sup>w</sup> (Eq. (31)); (2) peak frequency ωp; (3) energy* E *(Eq. (27)); and (4) approximation of the peak frequency evolution (Eq. (32)).*

On the whole, the evolutions of integral characteristics of the solution shown in **Figures 2** and **3** are smoother than those calculated by Chalikov [26]. It can be explained by multiple technical improvements of the numerical scheme and higher resolution.

The evolution of wave spectrum is shown in **Figure 4**.

The 2-D wave spectrum *S k*ð Þ , *l* 0≤*k*≤ *Mx*, �*My* ≤*l* ≤ *My* � � averaged over nine time intervals of length equal to Δ*t*≈ 500 was transferred to the polar coordinates *Sp*ð Þ� *ψ*,*r* ð Þ *π=*2≤ *ψ* ≤*π=*2, 0≤ *r*≤ *Mx* and then averaged over the angle *ψ* to obtain the 1-D spectrum*Sh*ð Þ*r* :

$$\mathcal{S}\_h(r) = \sum \mathcal{S}\_p(\boldsymbol{\varphi}, r) r \Delta \boldsymbol{\varphi}. \tag{33}$$

The angle *ψ* ¼ 0 coincides with the direction of wind *U*, Δ*ψ* ¼ *π=*180*:* Even the averaged over angle spectrum looks quite irregular and contains multiple holes and peaks. The spectra are smoothed.

The two-dimensional wave spectra are shown in **Figure 5,** where the log <sup>10</sup>ð Þ *S*ð Þ *ψ*,*r* averaged over the successive eight periods of length Δ*t* ¼ 500 is given.

The first panel with a mark 0 refers to the initial conditions. The pictures well characterize the downshifting and angle spreading of spectrum due to the nonlinear interactions.

As seen, each spectrum consists of separated peaks and holes. This phenomenon was observed and discussed by Chalikov et al. [7]. The same results were obtained by Chalikov [26]. The repeated calculations with different resolutions showed that such structure of 2-D spectrum is typical. The locations of peaks cannot be explained by the fixed combination of interacting modes, since in different runs (with the same initial conditions but a different set of phases for the modes), the peaks are located in different locations in the Fourier space.

### **Figure 4.**

*The wave spectra Sh*ð Þ*r integrated over angle ψ in the polar coordinates and averaged over the consequent intervals of length about 500 units of the nondimensional time* t *(thin curves). The spectra are growing and shifting from right to left. Thick dashed curve is the dependence S* � *<sup>ω</sup>*�<sup>4</sup>*; dotted curve corresponds to S* � *<sup>ω</sup>*�<sup>5</sup> *.*

### **Figure 5.**

*Sequence of 2-D images of* lоg10ð Þ *S r*ð Þ , *ψ averaged over the consequent eight periods of length* Δ*t* ¼ 500*. The numbers indicate the period of averaging (the first panel marked 0 refers to the initial conditions). The spectra are normalized by the maximum value of spectrum 8. The horizontal axis corresponds to the wave numbers r* ¼ j j *k ; the angles are shown by rays.*

It is interesting to note that while increasing resolution, the patches with low energy extend. It can be supposed that the current and higher resolutions are excessive, and the process can be simulated with a lower resolution. This statement *High-Resolution Numerical Simulation of Surface Wave Development under the Action of Wind DOI: http://dx.doi.org/10.5772/intechopen.92262*

may be too optimistic, but it can be supported by the following arguments. The multi-mode wave mechanics is different from the multi-scale turbulent motion. The modeling of turbulence at increase of resolution just allows reproducing more details of motion. The increase of resolution in a wave model introduces other wave modes with different phase velocities. Due to dispersion, the solutions (i.e., the evolution of surface) in these two cases will be completely different. It means that the solution does not converge with increase of resolution, which makes no sense.

The situation can be saved if upon reaching the optimal resolution, the new added positions for the modes will not obtain the energy and not participate in solution. The existence of such effect should be carefully validated with the exact wave model. If this effect does not exist, it means that the results of simulations depend completely on the resolution, the reliable simulation of individual evolution of wave field being, in fact, impossible.

The method of calculation of the simulated one-dimensional input and dissipation spectra was described by Chalikov [26]; still, it will be explained here once again, though briefly.

The evolution of the integrated over angle *ψ* wave spectrum *Sh*ð Þ*r* can be described with the equation:

$$\frac{d\mathbb{S}\_h(r)}{dt} = I(r) + D\_t(r) + D\_b(r) + N(r),\tag{34}$$

where *I r*ð Þ, *Dt*ð Þ*r* , *Db*ð Þ*r* and *N r*ð Þ are the spectra of the input energy, tail dissipation, breaking dissipation, and the rate of nonlinear interactions. All of the spectra shown below were obtained by transformation of the 2-D spectra into a polar coordinate ð Þ *ψ*,*r* and then integrated over the angles *ψ* within the interval ð Þ �*π=*2, *π=*2 . The spectra can be calculated using an algorithm similar to Eq. (29) for integral characteristics. For example, the spectrum of the energy input *I k*ð Þ , *l* is calculated as follows:

$$I(k,l) = \left(\mathbb{S}\_c^{t+\Delta t}(k,l) - \mathbb{S}\_c^t(k,l)\right) / \Delta t,\tag{35}$$

where *Sc kx*, *ky* � � is a spectrum of the columnar energy calculated by the relation:

$$S\_{\varepsilon}(k,l) = \frac{1}{2} \left( h\_{k,l}^2 + h\_{-k,-l}^2 + \int\_{-H}^0 (u\_{k,l}^2 + u\_{-k,-l}^2 + v\_{k,l}^2 + v\_{-k,-l}^2 + w\_{k,l}^2 + w\_{-k,-l}^2) d\zeta \right) \tag{36}$$

where the grid values of velocity components *u*, *v*, *w* are calculated by the relations:

$$
\mu = \varrho\_{\xi} + \varrho\_{\zeta}\eta\_{\xi}, \quad \upsilon = \varrho\_{\vartheta} + \varrho\_{\zeta}\eta\_{\vartheta}, \quad \ w = \varrho\_{\zeta}, \tag{37}
$$

and *uk*,*<sup>l</sup>*, *vk*,*<sup>l</sup>* and *wk*,*<sup>l</sup>* are the real Fourier coefficients, while for the negative indices—the imaginary ones.

For calculation of *I k*ð Þ , *l* , the fictitious time steps Δ*t* are made only with a term responsible for the energy input, that is, the surface pressure *p*. The spectrum *I k*ð Þ , *l* was averaged over the periods Δ*t*≈ 500, then transformed into a polar coordinate system and integrated in the Fourier space over the angles *ψ* within the interval ð Þ �*π=*2, *π=*2 *:* Such procedure was used for calculation of all the terms in the right side of Eq. (34). In the current version of the model, the calculations of integral (28) and spectral (34) transformations were combined with the calculations of the right sides of Eqs. (4) and (5).

The rates of transformation of spectrum are shown in **Figure 6**. The integral term describing the nonlinear interaction *N* in Eq. (28) is small (as compared with the local values of *Nk*,*l*), but the magnitude of spectrum *N r*ð Þ is comparable with the input *I r*ð Þ and dissipation *Dt*ð Þ*r* and *Db*ð Þ*r* terms (panel 1 in **Figure 5**). The shape of spectrum *N r*ð Þ confirms prediction of the quasi-linear theory [54, 55]. At the low wave number slope of the spectrum, the nonlinear influx of energy is positive, while at the opposite slope, it is negative. This process produces shifting of spectrum to the lower wave number (downshifting). The input of energy due to the nonlinear interactions is observed in a high frequency part of spectrum, which also agrees with Hasselmann's theory. Note that the nonlinear interactions also produce widening of spectrum.

The spectral distribution of the energy input from wind *I r*ð Þ (panel 2 in **Figure 6**) is in general similar to wave spectrum since it depends linearly on the spectral density (**Figure 3**). The dissipation rate *Db*ð Þ*r* is negative (panel 3), and its minimum is shifted a little to higher frequencies from the wave spectrum peak. The tail dissipation (Panel 4) is smaller by two orders than the other terms, but it plays an important role of supporting numerical stability.

*The rates of transformation spectra multiplied by* <sup>10</sup><sup>9</sup>*: (1) nonlinear interaction N r*ð Þ *(1); (2) input energy I r*ð Þ*; (3) breaking dissipation Db*ð Þ*<sup>r</sup> ; (4) tail dissipations Dt*ð Þ*<sup>r</sup> multiplied by* 1011*. All spectra are obtained by transformation of the 2-D spectra into the polar coordinate* ð Þ *ψ*,*r and then integrated over the angles ψ within the interval* ð Þ �*π=*2, *π=*2 *.*

*High-Resolution Numerical Simulation of Surface Wave Development under the Action of Wind DOI: http://dx.doi.org/10.5772/intechopen.92262*

### **Figure 7.**

*The spectral distribution of the residual input of energy dSh*ð Þ*r =dt in the energy-containing part of spectrum for different stages of wave development. The numbered vertical lines indicate positions of the spectral peak. The second set of numbers shows the peaks of the corresponding residual input.*

The residual rate of transformation of spectrum *dSh*ð Þ*r =dt* averaged over eight consequent periods is shown in **Figure 7**. The numbers in the top part of panel indicate the averaged wave number of the spectral peak. The second set of numbers refers to the corresponding spectrum of the residual input of energy. As seen, the maximum of the input energy is located to the left of the spectral peak, that is, on a low-wavenumber spectral slope. The obtained energy causes downshifting of spectrum and supports the shape of a high wavenumber slope and spectral tail. In the equilibrium regime, all the incoming energy are consumed for supporting the shape of the entire spectrum.

However, the dynamics of the tale is not adiabatic, that is, it is not completely controlled by the spectral energy cascade since the input of energy due to the

### **Figure 8.**

*The rates of transformation spectra multiplied by* 109 *for the last period: (1) nonlinear interaction N r*ð Þ *(1); (2) input energy; (3) breaking dissipation Db*ð Þ*<sup>r</sup> ; (4) tail dissipations* 1010*Dt*ð Þ*<sup>r</sup> ; and (5) balance of all the terms.*

nonlinear interaction competes with the energy input from wind (curve 2) and breaking dissipation (curves 3 and 4). The tail dissipation (curve 4) is small and concentrated in the vicinity of the cut wave number. The input and dissipation in the spectral tail are nearly in balance (curve 4) (**Figure 8**).
