# Evaluating instantaneous crest phase speed in multi-directional seas through Hilbert transform

Ying Wang, Guillaume Ducrozet

Nantes Université, École Centrale Nantes, CNRS, LHEEA, UMR 6598, France

Email: ying.wang@ec-nantes.fr, guillaume.ducrozet@ec-nantes.fr

## 1 Introduction

The definition of an appropriate wave-breaking onset criterion is essential for the development of highly nonlinear wave solvers. Barthelemy et al. (2018) proposed to use the ratio of the fluid velocity and crest phase speed, which has been demonstrated to be relevant in several configurations such as deep-water or depth-induced breaking. Then, accurately and efficiently computing the instantaneous wave crest phase speed is crucial. However, it appears challenging, particularly in multi-dimensional sea states. Seiffert et al. (2017) has demonstrated that the Hilbert Transform Method (HTM) efficiently captures the wave celerity for unidirectional irregular sea states. The present study extends the HTM to multi-directional sea states. A comparison between the HTM and the conventional Crest Tracking Method (CTM) reveals encouraging preliminary results. The numerical wave model adopted in the present work, HOS-NWT, is a computationally efficient open-source Numerical Wave Tank (NWT) based on the High-Order Spectral method, which has been extensively validated (Ducrozet et al. 2012).

## 2 Methodology

The Crest Tracking Method (CTM) tracks every wave crest at every time step. Denoting $(x_c(t), y_c(t))$ the crest locations, the corresponding phase velocity components in the x and y-direction $(c_x$ and $c_y$ respectively) are computed by taking the temporal derivative of those locations. The total wave crest phase velocity is therefore given by

$$c_n = \sqrt{(c_x)^2 + (c_y)^2} = \sqrt{(dx_c/dt)^2 + (dy_c/dt)^2} \quad (1)$$

Numerically, a least square method in time is utilized to evaluate the crest phase speed from the successive individual crest locations {$x_c(t_{i-n}), \dots, x_c(t_{i-1}), x_c(t_i)$}. Those are obtained by zero-crossing analysis together with a specific procedure to identify the wave associated with a given crest. Additionally, to reduce numerical oscillations, a time threshold $T_{th}$, is introduced to define the minimum time window for the least squares fitting. This defines the minimum duration of an individual wave before computing its crest phase speed. This is taken as $T_{th} = T_c/16$, with $T_c$ the characteristic period of the wavefield studied. The CTM is the most straightforward way to obtain the instantaneous crest phase speed. However, it suffers from its numerical complexity, which makes it difficult to use for large-scale problems in multi-directional seas.

The so-called Hilbert Transform Method (HTM) is using the partial Hilbert transform in space to compute the instantaneous wave number. Then, associated with the linear dispersion relation, it is possible to evaluate the instantaneous phase speed. Seiffert et al. (2017) has shown the efficiency of utilizing HTM to compute the crest phase speed
in unidirectional irregular sea states. However, in short-crested seas, the evaluation of the instantaneous phase function is more complicated. In the following, the temporal variable is treated as a constant since we only deal with the spatial variables, looking for information on the field at every time instant $t_i$. Therefore the three-dimensional wave elevation $\eta(x, y, t)$ reduces to two spatial dimensions $\eta(x, y, t_i)$. The partial and total spatial Hilbert transforms of the wave elevation $\eta(x, y, t_i)$ are defined as

$$ \mathbb{H}_x[\eta(x, y, t_i)] = \frac{1}{\pi} P \int_{-\infty}^{\infty} \frac{\eta(x', y, t_i)}{x - x'} dx' \quad ; \quad \mathbb{H}_y[\eta(x, y, t_i)] = \frac{1}{\pi} P \int_{-\infty}^{\infty} \frac{\eta(x, y', t_i)}{y - y'} dy' \quad (2) $$

$$ \mathbb{H}_{xy}[\eta(x, y, t_i)] = \frac{1}{\pi^2} P \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \frac{\eta(x', y', t_i)}{(x - x')(y - y')} dx' dy' \quad (3) $$

where $P$ is the Cauchy principal value of the integral, $\mathbb{H}_{xy}[\eta]$ denotes the total two-dimensional spatial Hilbert transform, $\mathbb{H}_x[\eta]$ and $\mathbb{H}_y[\eta]$ are the partial spatial Hilbert transforms with respect to $x$ and $y$, respectively. Omitting the time index $t_i$ for conciseness, for the two-dimensional signals, four complex signals are constructed:

$$ \Psi_1(x, y) = \eta(x, y) - \mathbb{H}_{xy}[\eta(x, y)] + j(\mathbb{H}_x[\eta(x, y)] + \mathbb{H}_y[\eta(x, y)]) \quad (4) $$

$$ \Psi_2(x, y) = \eta(x, y) + \mathbb{H}_{xy}[\eta(x, y)] - j(\mathbb{H}_x[\eta(x, y)] - \mathbb{H}_y[\eta(x, y)]) \quad (5) $$

$$ \Psi_3(x, y) = \eta(x, y) + \mathbb{H}_{xy}[\eta(x, y)] + j(\mathbb{H}_x[\eta(x, y)] - \mathbb{H}_y[\eta(x, y)]) \quad (6) $$

$$ \Psi_4(x, y) = \eta(x, y) - \mathbb{H}_{xy}[\eta(x, y)] - j(\mathbb{H}_x[\eta(x, y)] + \mathbb{H}_y[\eta(x, y)]) \quad (7) $$

where $j$ is the unit imaginary number, $j = \sqrt{-1}$. Eqs. (4~7) are rewritten as

$$ \Psi_1(x, y) = \Psi_4^*(x, y) = A_1(x, y) \exp\{j\varphi_1(x, y)\} \quad (8) $$

$$ \Psi_3(x, y) = \Psi_2^*(x, y) = A_2(x, y) \exp\{j\varphi_2(x, y)\} \quad (9) $$

where the superscript * represents the complex conjugate, $A_1(x, y)$ and $A_2(x, y)$ are the instantaneous amplitudes defined by

$$ A_1(x, y) = \sqrt{(\eta(x, y) - \mathbb{H}_{xy}[\eta(x, y)])^2 + (\mathbb{H}_x[\eta(x, y)] + \mathbb{H}_y[\eta(x, y)])^2} \quad (10) $$

$$ A_2(x, y) = \sqrt{(\eta(x, y) + \mathbb{H}_{xy}[\eta(x, y)])^2 + (\mathbb{H}_x[\eta(x, y)] - \mathbb{H}_y[\eta(x, y)])^2} \quad (11) $$

$\varphi_1(x, y)$ and $\varphi_2(x, y)$ are the instantaneous phases defined by

$$ \varphi_1(x, y) = \arctan \left( \frac{\mathbb{H}_x[\eta(x, y)] + \mathbb{H}_y[\eta(x, y)]}{\eta(x, y) - \mathbb{H}_{xy}[\eta(x, y)]} \right) \quad (12) $$

$$ \varphi_2(x, y) = \arctan \left( \frac{\mathbb{H}_x[\eta(x, y)] - \mathbb{H}_y[\eta(x, y)]}{\eta(x, y) + \mathbb{H}_{xy}[\eta(x, y)]} \right) \quad (13) $$

The free surface elevation can be written as a function of those amplitudes and phases

$$ \eta(x, y) = \frac{1}{2} \Re [A_1(x, y) \exp\{j\varphi_1(x, y)\} + A_2(x, y) \exp\{j\varphi_2(x, y)\}] \quad (14) $$

$$ \eta(x, y) = \frac{1}{2} A_{12}(x, y) \cos(\varphi_1(x, y) + \varphi_2(x, y) + \varphi_{12}(x, y)) \quad (15) $$

where $\Re[\cdot]$ is the function to take the real part. $A_{12}$ represents an amplitude modulation and $\varphi_{12}$ a new phase, which are respectively defined by

$$ A_{12}(x, y) = \sqrt{(A_1(x, y))^2 + (A_2(x, y))^2 + 2A_1(x, y)A_2(x, y) \cos(\varphi_1(x, y) - \varphi_2(x, y))} $$

$$ \varphi_{12}(x, y) = - \arctan \left( \frac{A_2 \sin(\varphi_1) + A_1 \sin(\varphi_2)}{A_2 \cos(\varphi_1) + A_1 \cos(\varphi_2)} \right) $$
Based on Eq.(15), the instantaneous phase function $\varphi(x, y)$ of the two-dimensional wave elevation $\eta(x, y)$ is

$$ \varphi(x, y) = \varphi_1(x, y) + \varphi_2(x, y) + \varphi_{12}(x, y) \quad (16) $$

This allows for the definition of the instantaneous wave number components in $x$- and $y$-directions ($k_x$ and $k_y$ respectively) using the corresponding partial derivatives

$$ k_x(x, y, t_i) = \frac{\partial \varphi(x, y, t_i)}{\partial x} \quad ; \quad k_y(x, y, t_i) = \frac{\partial \varphi(x, y, t_i)}{\partial y} \quad (17) $$

Therefore, the local wavenumber $k_n$ and direction $\theta_n$ can be evaluated as

$$ k_n(x, y, t_i) = \sqrt{k_x^2 + k_y^2} \quad ; \quad \theta_n(x, y, t_i) = \arctan \left( \frac{k_y(x, y, t_i)}{k_x(x, y, t_i)} \right) \quad (18) $$

Finally, it is possible to compute the instantaneous local crest phase speed using the linear dispersion relation, $h$ being the water depth, and $g$ the gravitational acceleration

$$ c_n(x, y, t_i) = \sqrt{\frac{g \cdot \tanh (k_n(x, y, t_i) \cdot h)}{k_n(x, y, t_i)}} \quad (19) $$

# 3 Results and Conclusion

A comparison between the CTM and the HTM is conducted with the simulation of a non-breaking directional focusing wave with HOS-NWT. We use the configuration reported in Brandini & Grilli (2002). The length, width, and water depth of the NWT are set as $L_x = 10.0$ m, $L_y = 16.0$ m and $h = 1.0$ m. The wavemaker motion is defined through the surface elevation at its location $x = 0$

$$ \eta(x = 0, y, t) = \sum_{l=1}^{N_{\theta}} a_l \cos \left[ k(y \sin \theta_l - x_f \cos \theta_l - y_f \sin \theta_l) - \omega t \right] \quad (20) $$

where $N_{\theta}$ is the number of directions, $a_l$, $k$, $\omega$ are wave amplitude, wave number, and angular frequency ($\omega = 2\pi f$). The focusing location is defined by $(x_f, y_f)$. Besides, $\theta_l$ is the propagating angle of each wave component and is uniformly distributed in the range $[-\theta_{max}, \theta_{max}]$. The number of wave components is set to $N_{\theta} = 4$, maximum angle $\theta_{max} = 25^{\circ}$, amplitude $a_l = 0.02$ m, wave number $k = 1.6866 \text{ m}^{-1}$, and focusing location $x_f = 7.45$ m, $y_f = 8.0$ m. In addition, the spatial discretization $N_x = 120(k_{xmax} \approx 22k)$, $N_y = 40$, and the HOS order is set as $M = 5$ in the HOS-NWT solver.

As an illustration, the wave crest phase velocity vectors, calculated by the CTM and the proposed HTM at $t = 6.80$ s, are presented in Fig.(1). The HTM computes a phase velocity at each spatial point, displayed here at the crests for each $y$-transect to validate the direction evaluation. Visually, at the crest location, the evaluation by the HTM shows encouraging results compared with the reference CTM. The detailed comparison of the crest phase speed at the maximum wave elevation among the whole surface is displayed in Fig.(1). The two methods overall agree well, with a far simpler implementation and computational effort for HTM. The discrepancies may be attributed to the nonlinearity, the HTM using linear dispersion relation. However, the computation by the CTM depends on a least square fitting, which may also lead to discrepancies. In contrast, the evaluation
<figure><img src="image_1.png" /></figure>

Figure 1: (Left) Illustration of the crest (green triangles) phase speed vector at $t = 6.8s$. Comparison between CTM (green arrows) and HTM (purple arrows); magenta lines represent the crest line. (Right) Comparison between $c_{ntrk}$, evaluated by the CTM, and $c_{nhil}$, computed by the HTM, at the location of the maximum wave elevation.

by HTM is not affected by the previous temporal information of the crest location and can be understood as an instantaneous methodology.

Overall, implementing the Hilbert transform to compute the crest phase speed is successfully extended to multi-directional sea states. The new methodology was tentatively validated with the simulation of a directional focusing wave. Compared to the CTM, the HTM is a computationally efficient, stable, and fully instantaneous method. In addition, the HTM can visualize the wave phase speed in the whole domain instead of only at the wave crest. As aforementioned, further validation is still required and will be presented at the workshop, including different types of directional wave conditions: frequency focusing as well as irregular sea states. This study is the preliminary step before the possible use of HTM in a wave-breaking model for nonlinear potential flow solvers.

## References

Barthelemy, X., Banner, M., Peirson, W., Fedele, F., Allis, M. & Dias, F. (2018), ‘On a unified breaking onset threshold for gravity waves in deep and intermediate depth water’, *Journal of Fluid Mechanics* **841**, 463–488.

Brandini, C. & Grilli, S. T. (2002), Three-dimensional wave focusing in fully nonlinear wave models, in ‘Ocean Wave Measurement and Analysis (2001)’, pp. 1102–1111.

Ducrozet, G., Bonnefoy, F., Le Touzé, D. & Ferrant, P. (2012), ‘A modified high-order spectral method for wavemaker modeling in a numerical wave tank’, *European Journal of Mechanics-B/Fluids* **34**, 19–34.

Seiffert, B. R., Ducrozet, G. & Bonnefoy, F. (2017), ‘Simulation of breaking waves using the high-order spectral method with laboratory experiments: Wave-breaking onset’, *Ocean Modelling* **119**, 94–104.