4. Theory

In this chapter we define the wave formulation as applied in the spectral_wave_data API.

The wave generator may apply other formulations but is required to do eventual output transformations to comply with the definitions described in this document. Such modifications are typical minor, limited to apply sign shifts, complex conjugate etc. in output statements.

The wave generator producing the wave field do not consider the location, orientation and speed of eventual structures. As a consequence, the produced ambient wave field can be reused for assessing different structural configurations.

4.1. The SWD coordinate system

The spectral_wave_data (SWD) wave field is described in an earth fixed right-handed Cartesian coordinate system \(\mathbf{x}=(x,y,z)\) where \(z=0\) coincides with the calm free surface and the \(z\)-axis is pointing upwards. This is the coordinate system applied in the SWD file.

For long crested seas, waves are assumed to propagate in the positive \(x\)-direction.

For short crested seas, waves may propagate in all directions. Consequently, the orientation of the \(x\)-direction is not specified by the API. However, if there is a time independent main propagation direction it is expected to be in the positive \(x\)-direction. In case of varying bathymetry this expectation may not apply.

4.2. The application defined wave coordinate system

The spectral_wave_data API is designed to be applied in a wide range of application programs. Consequently, an earth fixed, right-handed, Cartesian wave coordinate systems \(\mathbf{\bar{x}}=(\bar{x},\bar{y},\bar{z})\) related to the application program may differ from the SWD coordinate system. It is assumed that \(\bar{z}=0\) coincides with the calm free surface and the \(\bar{z}\)-axis is pointing upwards.

The relation between the SWD and the wave coordinate system applied in the application program is defined by the 3 spatial shift parameters \(x_0\), \(y_0\) and \(\beta\). The location of the application origin is \((x_0, y_0, 0)\) observed from the SWD system. The \(x\)-axis is rotated \(\beta\) relative to the \(\bar{x}\)-axis.

These relations are expressed as

\[x - x_0 = \bar{x}\cos\beta + \bar{y}\sin\beta, \qquad y - y_0 = -\bar{x}\sin\beta + \bar{y}\cos\beta\]
\[z = \bar{z}, \qquad t - t_0 = \bar{t}, \qquad t_0 \ge 0\]

The shift parameter \(t_0\) relates the time \(t\) defined in the SWD system relative to the time \(\bar{t}\) in the application program.

_images/corsys.png

The relation between the SWD \(\mathbf{x}\) and the earth fixed application program coordinate system \(\mathbf{\bar{x}}\). In this example \(x_0\), \(y_0\) and \(\beta\) are all positive quantities.

4.3. The spectral formulation

Assuming an ocean of ideal fluid, the incident velocity potential \(\phi(\mathbf{x}, t)\) and the corresponding single valued surface elevation \(\zeta(x, y, t)\) at time \(t\) are expressed with respect to the SWD coordinate system as

\[\phi(\mathbf{x}, t)= \sum_{j=1}^n \mathcal{Re} \Bigl\{c_j(t)\, \psi_j(\mathbf{x})\Bigr\}, \qquad \zeta(x, y, t)= \sum_{j=1}^n \mathcal{Re} \Bigl\{h_j(t)\, \xi_j(x, y)\Bigr\}\]

where both sets of \(n\) explicit shape functions \(\psi_j(\mathbf{x})\) and \(\xi_j(x, y)\) are orthogonal on the global considered space. The functions \(\psi_j(\mathbf{x})\) are harmonic. The number of shape functions \(n\) defines a time independent spatial resolution of the wave field. Consequently, constrained by the given resolution \(n\) a wide range of wave fields can be described given the shape functions and the temporal amplitudes \(c_j(t)\) and \(h_j(t)\). Obvious exceptions are overturning waves and viscous flows, but a large class of important nonlinear wave trains may be described in addition to all the classical perturbation wave models (Airy, Stokes, Stream waves, etc.).

The shape functions and temporal amplitudes may in general be complex functions. We let \(\mathcal{Re}\{\alpha\}\) and \(\mathcal{Im}\{\alpha\}\) denote the real and imaginary part of a complex number \(\alpha\).

4.4. Spectral kinematics

Given the wave potential as defined above, the API may evaluate kinematics as described in the following table.

Symbol

API

Note

\(\phi\)

phi()

The velocity potential is often applied as Dirichlet conditions in
potential flow solvers.

\(\varphi\)

stream()

The stream function is orthogonal to the velocity potential. \(\frac{\partial \phi}{\partial x}=\frac{\partial\varphi}{\partial z}\)
and \(\frac{\partial\phi}{\partial z} = -\frac{\partial\varphi}{\partial x}\) The stream function is only relevant for long crested
waves. For short crested seas we define \(\varphi\equiv 0\) for all \(\mathbf{x}\) and \(t\).

\(\frac{\partial\phi}{\partial t}\)

phi_t()

is the time rate of change of \(\phi\) at an earth fix location. It is
sometimes applied as Dirichlet conditions in potential flow solvers. This
term is an important part of the pressure field. For linear models it is
proportional to the dynamic pressure.

\(\nabla\phi\)

grad_phi()

The fluid particle velocity

\(\nabla\nabla\phi\)

grad_phi_2nd()

The second order spatial gradients of \(\phi\) are sometimes applied as
boundary conditions in potential flow solvers. They are also needed for
calculating particle accelerations.

\(\frac{\partial\nabla\phi}{\partial t}\)

acc_euler()

The Euler acceleration. (Earth fixed location)

\(\frac{d\nabla\phi}{dt}\)

acc_particle()

The fluid particle acceleration

\(p\)

pressure()

The fluid pressure including all terms in Bernoulli’s equation.

\(\zeta\)

elev()

The surface elevation

\(\frac{\partial\zeta}{\partial t}\)

elev_t()

The time rate of change of surface elevation (At fixed horizontal location)

\(\nabla\zeta\)

grad_elev()

Spatial gradients of the surface elevation is often applied in boundary
conditions

\(\nabla\nabla\zeta\)

grad_elev_2nd()

The second order spatial gradients of the surface elevation are sometimes
applied in boundary conditions for potential flow solvers.

4.5. Shape classes

We define a shape class as a specific set of parameterized shape functions \(\psi_j(\mathbf{x})\) and \(\xi_j(x, y)\) as explained above. A wave field is assumed to be described by a single shape class.

Because \(\psi_j(\mathbf{x})\) is harmonic we may in general write

\[\psi_j(\mathbf{x})= X_j(x)\, Y_j(y)\, Z_j(z)\]

In the current version, the official API implements the following shape classes

  • Shape 1 For long crested waves propagating in infinite water depth.

  • Shape 2 For long crested waves propagating in constant finite water depth.

  • Shape 3 For long crested waves propagating in infinite, constant or varying water depth.

  • Shape 4 For short crested waves propagating in infinite water depth.

  • Shape 5 For short crested waves propagating in constant finite water depth.

  • Shape 6 For a general set of Airy waves. Infinite or constant water depth.

4.6. Temporal amplitudes

Applying a shape class as defined above, the kinematics at time \(t\) is well defined given the set of associated temporal amplitudes \(c_j(t)\) and \(h_j(t)\). The wave generators are required by this API to provide discrete values of these functions at equidistant time intervals. In order to obtain consistent smooth kinematics in application programs, it is also required that \(\frac{dc_j(t)}{dt}\) and \(\frac{dh_j(t)}{dt}\) are provided by the wave generators at the same time intervals. In order to minimize errors in application programs the latter functions should be provided directly from the kinematic and dynamic free surface conditions and not by using temporal finite difference schemes:

\[\frac{\partial\phi}{\partial t} = \sum_{j=1}^n \mathcal{Re} \Bigl\{\frac{dc_j(t)}{dt}\, \psi_j(\mathbf{x})\Bigr\} = f_d(\phi, \nabla\phi, \ldots)\]
\[\frac{\partial\zeta}{\partial t} = \sum_{j=1}^n \mathcal{Re} \Bigl\{\frac{dh_j(t)}{dt}\, \xi_j(x, y)\Bigr\} = f_k(\phi, \nabla\phi, \ldots)\]

where the functions \(f_d()\) and \(f_k()\) follows from the applied dynamic and kinematic free surface conditions and are calculated using spatial gradients only.

Note

  • For spectral solvers \(\frac{dc_j(t)}{dt}\) and \(\frac{dh_j(t)}{dt}\) are obtained from the right-hand-side of the temporal ordinary-differential-equation system.

  • Only a few time steps of \(c_j(t)\), \(h_j(t)\), \(\frac{dc_j(t)}{dt}\) and \(\frac{dh_j(t)}{dt}\) needs to be in memory for evaluating the kinematics at any time instance.

  • The presence of explicit values of \(\frac{dc_j(t)}{dt}\) and \(\frac{dh_j(t)}{dt}\) ensures that the kinematics produced by the API comply with the physical boundary conditions at the specified instances.

  • When applying interpolation schemes including constraints of derivatives, the issue of over-shooting is reduced.

4.6.1. Consistent interpolation

Given the value and slope of \(c_j(t)\) and \(h_j(t)\) at equidistant time instances \(t_{i}, t_{i+1},\ldots\) it is important that the predicted kinematics between these intervals are consistent and sufficiently smooth for proper application in other solvers. E.g. if the potential \(\phi\) is assumed linear between two time steps the pressure distribution will be discontinuous at every new time interval due to the pressure component \(\frac{\partial\phi}{\partial t}\). This is not acceptable.

In order to deal with simulations over a long time range it is important that the constant time intervals \(\Delta t= t_{i+1} - t_i\) can be relatively large. Even slightly larger than in the program providing the kinematics. The program that applies the kinematics may have order of magnitude smaller time steps like in slamming CFD simulations. The kinematics provided by this API should not introduce numerical noise in such calculations.

Let \(f_i\) and \(\frac{df_i}{dt}\) denote \(c_j(t_i)\) and \(\frac{dc_j(t_i)}{dt}\) or \(h_j(t_i)\) and \(\frac{dh_j(t_i)}{dt}\) respectively. When the application program updates the actual time \(t\) to an arbitrary value, the API will once approximate the values of \(f(t)\) and \(\frac{df(t)}{dt}\) using one of the interpolation schemes described below. Then any kinematics at time \(t\) can be evaluated using these current spectral amplitudes.

Two alternative interpolation schemes are provided.

4.6.1.1. \(C^1\) continuous spline

Given the spectral data at two subsequent time steps \(t_i\) and \(t_{i+1}\)

\[f_i, f_{i+1}, \frac{df_i}{dt}, \frac{df_{i+1}}{dt}\]

we obtain a cubic spline which is \(C^1\) continuous at the boundary of each time interval.

\[f(t) = (1-\delta)f_i + \delta f_{i+1} + \delta(1-\delta)(a(1-\delta) + b \delta), \qquad t\in[t_i, t_{i+1}]\]
\[\frac{df}{dt}(t) = \frac{f_{i+1} - f_i}{\Delta t} + (1-2\delta)\frac{a(1-\delta) + b\delta}{\Delta t} + \delta(1-\delta)\frac{b-a}{\Delta t}\]
\[\delta = \frac{t-t_i}{\Delta t}, \qquad a = \Delta t\frac{df_i}{dt} - (f_{i+1} - f_i)\]
\[b = -\Delta t\frac{df_{i+1}}{dt} + (f_{i+1} - f_i)\]

This spline will produce a pressure field which is \(C^0\) continuous at the transition of each time step.

4.6.1.2. \(C^2\) continuous spline

Note

This is the default scheme in the API.

Some boundary value formulations require continuity of quantities like \(\frac{\partial^2\phi}{\partial t^2}\) and \(\frac{\partial^2\zeta}{\partial t^2}\). Since we know the first order derivatives from the dynamic and kinematic free surface conditions we construct a spline which is second order continuous at the transitions of time steps by first approximating the second order gradients at these points.

The following approximations

\[\frac{d^2 f_i}{dt^2} \approx 2\frac{f_{i-1} - 2 f_i + f_{i+1}}{\Delta t^2} - \frac{\frac{d f_{i+1}}{dt} - \frac{d f_{i-1}}{dt}}{2\Delta t}\]
\[\frac{d^2 f_{i+1}}{dt^2} \approx 2\frac{f_{i} - 2 f_{i+1} + f_{i+2}}{\Delta t^2} - \frac{\frac{d f_{i+2}}{dt} - \frac{d f_{i}}{dt}}{2\Delta t}\]

are obtained applying the central finite difference scheme as explained in this figure.

_images/interpolation5th.png

The green curve is the lowest order polynomial matching the values and slopes at \(t_{i-1}\), \(t_i\), and \(t_{i+1}\). The magenta curve is the lowest order polynomial matching the values and slopes at \(t_i\), \(t_{i+1}\), and \(t_{i+2}\). The 5th order polynomial \(f(t)\) matches the values and slopes at \(t_i\), and \(t_{i+1}\). In addition \(f(t)\) matches the curvature of the green curve at \(t_i\) and the curvature of the magenta curve at \(t_{i+1}\).

Knowing the value, slope and approximate curvature at \(t_i\) and \(t_{i+1}\) the lowest order unique spline \(f(t)\), and its derivative \(df/dt(t)\), on this interval follows.

\[f(t) = f_i +\sum_{i=1}^5 q_i\delta^i, \qquad \frac{d f}{dt}(t) = \frac{d f_i}{dt} + \frac{1}{\Delta t}\sum_{i=2}^5 i q_i\delta^{i-1} \qquad t\in[t_i, t_{i+1}]\]
\[\delta = \frac{t-t_i}{\Delta t} \in [0,1), \quad q_1 = \frac{d f_i}{dt} \Delta t, \quad q_2 = f_{i-1} - 2f_i + f_{i+1} + \Bigl(\frac{d f_{i-1}}{dt} - \frac{d f_{i+1}}{dt}\Bigr)\frac{\Delta t}{4}\]
\[q_3 = -3f_{i-1} - 3f_i + 5f_{i+1} +f_{i+2} - \Bigl( 3\frac{d f_{i-1}}{dt} + 23\frac{d f_i}{dt} + 13\frac{d f_{i+1}}{dt} + \frac{d f_{i+2}}{dt}\Bigr)\frac{\Delta t}{4}\]
\[q_4 = 3f_{i-1} + 7 f_{i} - 8f_{i+1}- 2f_{i+2} + \Bigl( 3\frac{d f_{i-1}}{dt} + 30\frac{d f_i}{dt} + 25\frac{d f_{i+1}}{dt} + 2\frac{d f_{i+2}}{dt}\Bigr)\frac{\Delta t}{4}\]
\[q_5 = -f_{i-1} - 3f_{i} + 3f_{i+1} + f_{i+2} - \Bigl( \frac{d f_{i-1}}{dt} + 11\frac{d f_i}{dt} + 11\frac{d f_{i+1}}{dt} + \frac{d f_{i+2}}{dt}\Bigr)\frac{\Delta t}{4}\]

At the first time step (\(i=1\)) of the simulation we apply the following “constant acceleration” padding of data when applying the formulas above.

\[f_0 = f_1 + \Bigl(\frac{d f_2}{dt} - 3\frac{d f_1}{dt}\Bigr)\frac{\Delta t}{2}\]
\[\frac{d f_0}{dt} = 2\frac{d f_1}{dt} - \frac{d f_2}{dt}\]

Similar we apply the following padding at the very last time step (\(i=n\)).

\[f_{n+1} = f_n - \Bigl(\frac{d f_{n-1}}{dt} - 3\frac{d f_n}{dt}\Bigr)\frac{\Delta t}{2}\]
\[\frac{d f_{n+1}}{dt} = 2\frac{d f_n}{dt} - \frac{d f_{n-1}}{dt}\]

It should be noted that the spline \(f(t)\) is constructed using only 4 time steps in memory. However, the curve is still \(C^2\) continuous at the transition of each time step.