**2. Equations**

The *nonstationary surface-following nonorthogonal* coordinate system is used as follows:

$$
\xi = \varkappa, \ \vartheta = \mathcal{y}, \ \zeta = \mathfrak{z} - \eta(\xi, \mathfrak{z}, \mathfrak{z}), \ \mathfrak{z} = \mathfrak{t}, \tag{1}
$$

where *η*ð Þ¼ *x*, *y*, *t η ξ*ð Þ , *ϑ*, *τ* is a moving periodic wave surface given by the Fourier series

$$\eta(\xi,\theta,\tau) = \sum\_{-M\_{\mathcal{X}} < k < M\_{\mathcal{X}} - M\_{\mathcal{I}} < l < M\_{\mathcal{I}}} \sum\_{l < M\_{\mathcal{I}}} h\_{k,l}(\tau) \Theta\_{k,l},\tag{2}$$

where *k* and *l* are the components of the wave number vector k; *hk*,*<sup>l</sup>*ð Þ*τ* are Fourier amplitudes for elevations *η ξ*ð Þ , *ϑ*, *τ* ; *Mx* and *My* are the numbers of modes in the directions *ξ* and *ϑ*, respectively, while Θ*<sup>k</sup>*,*<sup>l</sup>* are the Fourier expansion basis functions represented as the matrix:

$$\Theta\_{kl} = \begin{cases} \cos\left(k\xi + l\vartheta\right) & -M\_{\ge} \le k \le M\_{\ge}, \quad -M\_{\ge} < l < 0 \\\\ \cos\left(k\xi\right) & -M\_{\ge} \le k \le 0, \quad l = 0 \\\\ \sin\left(k\xi\right) & 0 \le k \le M\_{\ge}, \quad l = 0 \\\\ \sin\left(k\xi + l\vartheta\right) & -M\_{\ge} \le k \le M\_{\ge}, \quad 0 < l \le M\_{\ge} \end{cases} \tag{3}$$

The 3-D equations of potential waves in the system of coordinates (1) at *ζ* ≤0 take the following form:

$$
\eta\_{\tau} = -\eta\_{\xi}\rho\_{\xi} - \eta\_{\vartheta}\rho\_{\vartheta} + \left(\mathbf{1} + \eta\_{\xi}^{2} + \eta\_{\vartheta}^{2}\right)\Phi\_{\xi},\tag{4}
$$

$$\rho\_{\tau} = -\frac{1}{2} \left( \rho\_{\xi}^{2} + \rho\_{\theta}^{2} - \left( \mathbf{1} + \eta\_{\xi}^{2} + \eta\_{\theta}^{2} \right) \Phi\_{\zeta}^{2} \right) - \eta - p,\tag{5}$$

$$
\Phi\_{\xi\xi} + \Phi\_{\theta\theta} + \Phi\_{\zeta\zeta} = \Upsilon(\Phi),
\tag{6}
$$

where ϒ is the operator:

$$\Upsilon(\cdot) = 2\eta\_{\xi}(\cdot)\_{\xi\zeta} + 2\eta\_{\theta}(\cdot)\_{\theta\zeta} + \left(\eta\_{\xi\xi} + \eta\_{\theta\theta}\right)(\cdot)\_{\zeta} - \left(\eta\_{\xi}^2 + \eta\_{\theta}^2\right)(\cdot)\_{\zeta\zeta},\tag{7}$$

capital fonts Φ are used for the domain *ζ* <0, while the lower case *φ* refers to *ζ* ¼ 0*:* The term *p* in Eq. (5) describes the pressure on the surface *ζ* ¼ 0.

It is suggested by Chalikov et al. [7] that it is convenient to represent the velocity potential *<sup>φ</sup>* as a sum of the two components: a linear component <sup>Φ</sup>, *<sup>φ</sup>* <sup>¼</sup> <sup>Φ</sup>ð Þ *<sup>ξ</sup>*, *<sup>ϑ</sup>*, 0 � � and an arbitrary nonlinear component <sup>Φ</sup><sup>~</sup> , *<sup>φ</sup>*<sup>~</sup> <sup>¼</sup> <sup>Φ</sup><sup>~</sup> ð Þ *<sup>ξ</sup>*, *<sup>ϑ</sup>*, 0 � �:

$$
\rho = \overline{\rho} + \tilde{\rho}, \quad \Phi = \overline{\Phi} + \tilde{\Phi}. \tag{8}
$$

The linear component Φ satisfies Laplace equation:

$$
\overline{\Phi}\_{\xi\xi} + \overline{\Phi}\_{\theta\theta} + \overline{\Phi}\_{\zeta\zeta} = 0,\tag{9}
$$

with a known solution:

$$\overline{\Phi}(\xi,\theta,\zeta,\mathfrak{r}) = \sum\_{k,l} \overline{\rho}\_{k,l}(\mathfrak{r}) \exp\left(|k|\zeta\right) \Theta\_{k,l},\tag{10}$$

where j j *<sup>k</sup>* <sup>¼</sup> *<sup>k</sup>*<sup>2</sup> <sup>þ</sup> *<sup>l</sup>* <sup>2</sup> � �<sup>1</sup>*=*<sup>2</sup> , *φ<sup>k</sup>*,*<sup>l</sup>* are the Fourier coefficients of the surface linear potential *φ* at *ζ* ¼ 0Þ*:* The solution satisfies the following boundary conditions:

$$\begin{aligned} \mathfrak{\zeta} &= \mathbf{0} : & \overline{\Phi} &= \overline{\rho} \\ \mathfrak{\zeta} & \to -\infty : & \tilde{\Phi}\_{\zeta} & \to \mathbf{0} \end{aligned} \tag{11}$$

The nonlinear component satisfies the equation:

$$
\tilde{\Phi}\_{\xi\xi} + \tilde{\Phi}\_{\theta\theta} + \tilde{\Phi}\_{\zeta\zeta} = \Upsilon(\tilde{\Phi}) + \Upsilon(\overline{\Phi}), \tag{12}
$$

Eq. (12) is solved with the boundary conditions:

$$\begin{aligned} \mathfrak{c} = \mathbf{0} &: & \tilde{\Phi} = \mathbf{0} \\ \mathfrak{c} \to -\infty &: & \tilde{\Phi}\_{\zeta} \to \mathbf{0} \end{aligned} \tag{13}$$

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

The presentation (8) is not used for solution of the evolutionary Eqs. (4) and (5) because it does not provide any improvements of accuracy and speed.

Eqs. (4)–(6) are written in a nondimensional form by using the following scales: length *L,* where 2π*L* is a (dimensional) period in the horizontal direction; time *L*1/2*g* �1/2; and velocity potential *L*3*=*<sup>2</sup> *g*1*=*<sup>2</sup> (*g* is the acceleration of gravity). The pressure is normalized by the water density, so that the pressure scale is *Lg*. Eqs. (4)–(6) are self-similar to the transformation with respect to *L.* The dimensional size of the domain is 2*πL*, so the scaled size is 2*π:* All of the results presented in this paper are nondimensional. Note that the number of the Fourier modes can be different in the *x* and *y* directions. In this case, it is assumed that the two-length scales *Lx* and *Ly* are used. The nondimensional length of the domain in the *y*-direction remains equal to 2*π*, and a factor *r* ¼ *Lx=Ly* is introduced into the definition of a differential operator in the Fourier space.

The derivatives of a linear component Φ in (7) are calculated analytically. The scheme combines a 2-D Fourier transform method in the 'horizontal surfaces' and a second-order finite-difference approximation on the stretched staggered grid defined by the relation Δ*ζ <sup>j</sup>*þ<sup>1</sup> ¼ *χ*Δ*ζ <sup>j</sup>* (Δ*ζ* is a vertical step, while *j* ¼ 1 at the surface). The stretched grid provides an increase in accuracy of approximation for the exponentially decaying modes. The values of the stretching coefficient *χ* lie for different settings within the interval 1.01–1.20; in the current work, the value*γ* ¼ 1*:*2 was used at the number of levels *Lw* ¼ 10. Such poor resolution was possible to use because of the separation of the potential into a large linear and a small nonlinear part, so Eq. (12) was used only for calculation of a small correction for the potential. A high value of the stretching coefficient provided high resolution in the vicinity of surface for accurate calculations of the surface vertical derivative for the potential. The finite-difference second-order approximation of vertical operators in Eq. (12) on a nonuniform vertical grid is quite straightforward (see Ref. [8]). Eq. (12) is solved as Poisson equations with the iterations over the right-hand side by TDMA method [27]. At each time step, the iterations start with a right-hand side calculated at the previous time step. A relative accuracy of the solution in terms of the vertical derivative of the potential on the surface was equal to 10�6. The typical number of iterations was 2–5.

The accuracy of the adiabatic version of equations was validated by reproducing a moving Stokes wave with the steepness *AK* ¼ 0*:*40 (*А* is a half of trough-to-crest wave height; *K* ¼ 1 is the wave number of the first mode). An algorithm for calculation of Stokes wave with the prescribed accuracy was suggested by Chalikov and Sheinin [28]. The scheme based on the conformal coordinates is very effective: the calculations were carried out in 100 ms at notebook (2.10 GHz). The dependence of Stokes wave spectra on the wave number is shown in **Figure 1**. In fact, about 2000 curves obtained in the course of calculations, with the interval Δ*t* ¼ 1, were plotted. Due to improvement of the numerical scheme, the accuracy of reproduction of Stokes wave is considerably higher than for the scheme used in Refs. [7, 8].

As seen, up to *<sup>S</sup>* � <sup>10</sup>�<sup>12</sup> ð Þ *<sup>k</sup>* <sup>¼</sup> <sup>22</sup> , the spectra of Stokes waves remain with high accuracy the same as it was assigned in initial conditions. At higher frequencies, the random disturbances appear. Note that this validation is not trivial: even small inaccuracies in the numerical scheme cause a fast distortion of the spectra, like in the bottom part of **Figure 1**. We consider these results as a serious evidence of high accuracy of the adiabatic version of the model. The previous version of 3-D model [26] allowed carrying out a long simulation of Stokes wave not steeper than *AK* ¼ 0*:*30. The right-hand sides of Eqs. (4) and (5) were calculated with a use of Fourier transform method: the nonlinear terms were calculated at the extended grid with size 4*Mx* � 4*My* , and then by the inverse Fourier transform, they were returned

**Figure 1.** *The 2000 spectra of Stokes wave with steepness AK* ¼ 0*:*40 *as a function of wave number k superimposed on each other. The gray area in the right left corner is a computational noise.*

to the Fourier grid. The fourth-order Runge-Kutta scheme was used for integration in time. The equation for the potential was solved at each of the four substeps of time step.

The simulations described by Chalikov [26] were a first attempt to reproduce the development of wave field assigned in the initial conditions as a group of small waves at high wave number under the action of strong wind. The initial elevation was generated as a superposition of linear waves corresponding to JONSWAP spectrum [29] with random phases. The initial Fourier amplitudes for the surface potential were calculated by the formulas of the linear wave theory. The details of the initial conditions are of no importance because the initial energy level is quite low. The wave peak was placed to the wave number equal to 100. The wind velocity was assigned equal to4*c*100, where *c*<sup>100</sup> is a phase velocity of the 100th mode. A detailed description of the scheme and its validation is given in Refs. [7, 8].

The simulation described in the current paper was performed with a doubled resolution in both directions, with the improved numerical scheme for Poisson equations and modified parameters in the scheme for calculations of energy transitions.
