2.3.2. Numerical models (Part I)

2.3.2.1. Finite element formulation

This section presents the formulation based on the finite element method (FEM) to solve the system of the wave diffraction-radiation problem. This formulation has been developed to be used in conjunction with unstructured meshes. The use of unstructured meshes enhances geometry flexibility and speed ups the initial modelling time. Let \(Q_h*\) be the finite element space to interpolate functions, constructed in the usual manner. From this space, we can construct the subspace \(Q_{h\Phi}\), that incorporates the Dirichlet conditions for the potential \(\Phi\). The space of test functions, denoted by \(Q_h\), is constructed as \(Q_{h\Phi}\), but with functions vanishing on the Dirichlet boundary. The weak form of the problem can be written as follows:

Find \([\phi_h] \in Q_{h,\phi}\), by solving the discrete variational problem:

(2.434)\[\begin{split}\begin{aligned} & \int_{\Omega} \nabla v_h \cdot \nabla \phi_h d\Omega = \\ & \int_{\Gamma^B} v_h \cdot \widehat{\phi}_n^{\small{B}} d\Gamma + \int_{\Gamma^R} v_h \cdot \widehat{\phi}_n^{\small{R}} d\Gamma + \int_{\Gamma^{Z_0}} v_h \cdot \widehat{\phi}_n^{\small{Z_0}} d\Gamma + \int_{\Gamma^{Z_{-H}}} v_h \cdot \widehat{\phi}_n^{\small{Z_{-H}}} d\Gamma \\ \\ & \forall v_h \in Q_h \end{aligned}\end{split}\]

Where \(\widehat{\phi}_n^{\small{B}}\), \(\widehat{\phi}_n^{\small{R}}\), \(\widehat{\phi}_n^{\small{Z_0}}\) and \(\widehat{\phi}_n^{\small{Z_{-H}}}\) are the potential normal gradients corresponding to the Neumann boundary conditions on bodies, radiation boundary, free surface and bottom, respectively. At this point, it is useful to introduce the associated matrix form:

(2.435)\[\begin{aligned} & \overline{\overline{\boldsymbol{L}}} \boldsymbol{\Phi} = \boldsymbol{b}^B + \boldsymbol{b}^R + \boldsymbol{b}^{Z_0} + \boldsymbol{b}^{Z_{-H}} \end{aligned}\]

where \(\overline{\overline{\boldsymbol{L}}}\) is the standard laplacian matrix, and \(\boldsymbol{b}^B\), \(\boldsymbol{b}^R\), \(\boldsymbol{b}^{Z_0}\), and \(\boldsymbol{b}^{Z_{-H}}\) are the vectors resulting of integrating the corresponding boundary condition terms. Regarding the bottom boundary for the refracted and radiated potential, it is imposed naturally in FEM by \(\boldsymbol{b}^{Z_{-H}} = 0\). A detailed description of the FEM model implemented in SeaFEM can be found in [1]

2.3.2.2. Free surface boundary condition

Combining the kinematic and dynamic free surface boundary conditions, the free surface condition reads:

(2.436)\[\begin{aligned} & \frac{\partial^2 \phi}{\partial t^2} + g \frac{\partial\phi}{\partial z} + \frac{\partial}{\partial t} \left(\frac{P_{fs}}{\rho }\right) + \{ Q^1 \} = 0 \end{aligned}\]

and is implemented as a Neumann boundary condition that fulfils the flux boundary integral:

(2.437)\[\begin{aligned} & \boldsymbol{b}^{Z_0} = \overline{\overline{\boldsymbol{M}}}_{\Gamma^{Z_0}} \end{aligned}\]

where \(\overline{\overline{\boldsymbol{M}}}_{\Gamma^{Z_0}}\) is the corresponding boundary mass and \(Q^1\) are the transfer terms from the first-order to the second-order problem (see Eqs. (1-54)-(1-55)):

(2.438)\[\begin{aligned} & Q^1=\partial_t R^1 - S^1 \end{aligned}\]

(2.439)\[\begin{aligned} & R^1 = \eta^1 \frac{\partial}{\partial z} \left(\frac{\partial\phi^1}{\partial t}\right) + \zeta^1 \frac{\partial}{\partial z} \left(\frac{\partial\phi^1}{\partial t}\right) + \eta^1 \frac{\partial}{\partial z} \left(\frac{\partial\psi ^1}{\partial t}\right) + \frac{1}{2} \nabla \phi^1 \cdot \nabla \phi^1 + \nabla \psi ^1 \cdot \nabla \phi^1 \end{aligned}\]

(2.440)\[\begin{split}\begin{aligned} & S^1 = \frac{\partial\phi^1}{\partial x} \frac{\partial\eta^1}{\partial x} + \frac{\partial\phi^1}{\partial y} \frac{\partial\eta^1}{\partial y} + \frac{\partial\phi^1}{\partial x} \frac{\partial \zeta^1}{\partial x} + \frac{\partial\phi^1}{\partial y} \frac{\partial \zeta^1}{\partial y} + \frac{\partial \psi ^1}{\partial x} \frac{\partial\eta^1}{\partial x} + \frac{\partial \psi ^1}{\partial y} \frac{\partial\eta^1}{\partial y} \\ & -\eta^1 \frac{\partial^2 \phi^1}{\partial z^2} -\eta^1 \frac{\partial^2 \psi ^1}{\partial z^2} -\zeta^1 \frac{\partial^2 \phi^1}{\partial z^2} \end{aligned}\end{split}\]

where

(2.441)\[\begin{aligned} & (\partial_t R^1)^{n+1} = 25/12 (R^1 )^{n+1} - 4 (R^1 )^n + 3 (R^1)^{n-1} - 4/3 (R^1 )^{n-2} + 1/4 (R^1)^{n-3} \end{aligned}\]

Eq. (2-3) is discretized in time using the following numerical scheme:

(2.442)\[\begin{split}\begin{aligned} & \frac{ \phi^{n+1} -2\phi^n + \phi^{n-1} }{ \Delta t^2 } = -g\phi_z^n -\frac{1}{12\Delta t^2} g \left(\phi_z^{n+1} + 10 \phi_z^n + \phi_z^{n-1} \right) - \\ & \frac{ P_{fs}^{n+1} - P_{fs}^{n-1} }{ \rho 2\Delta t } + \Big\{ -\frac{1}{12\Delta t^2} \left( (Q^1)^{n+1} + 10(Q^1)^n + (Q^1)^{n-1} \right) \Big\} \end{aligned}\end{split}\]

where for the specific case where \(P_{fs} =0\), the above scheme becomes a fourth order compact Padé scheme.The free surface elevation is discretized in time using the following fourth order in time numerical scheme:

(2.443)\[\begin{split}\begin{aligned} & \eta^{n+1} = - \frac{1}{g\Delta t} \left( \frac{25}{12} \phi^{n+1} -4\phi^n + 3\phi^{n-1} - \frac{4}{3} \phi^{n-2} + \frac{1}{4} \phi^{n-3} \right) - \\ & - \frac{P_{fs}^{n+1}}{\rho g} \Big\{-(S^1)^{n+1} \Big\} \end{aligned}\end{split}\]

2.3.2.3. Radiation condition and wave absorption

Waves represented by \(\pi\) are born at the bodies and propagate in all directions away from the bodies. These waves have to either be dissipated or to be let go out the domain so they will not come back and interact with the bodies. Because of this, a Sommerfeld radiation condition at the edge of the computational domain is introduced:

(2.444)\[\begin{aligned} & \partial _t \phi + c\nabla \phi \cdot n_R = 0 & \text{in } \Gamma_R \end{aligned}\]

where \(\Gamma_R\) is the surface limiting the domain in the horizontal directions, and \(c\) is a prescribed wave phase velocity. This radiation condition will let waves moving at velocity \(c\) to escape out the domain. The numerical scheme used to implement the radiation condition is

(2.445)\[\begin{aligned} & \left( \phi_n^R \right)^{n+1} = -\frac{\phi^{n-1}-\phi^n}{c\Delta t} & \text{in } \Gamma_R \end{aligned}\]

However, waves with very different velocities will not be leaving the domain. Hence, wave absorption is also introduced into the dynamic free surface boundary condition by varying the pressure such that:

(2.446)\[\begin{aligned} & P_{fs} (\boldsymbol{x},t) = \kappa (\boldsymbol{x})\rho \frac{\partial \phi }{\partial z} \end{aligned}\]

Eq. (2-12) increases pressure when the free surface is moving upwards, while decreases the pressure when the free surface is moving downwards. By doing this, energy is transferred from the waves to the atmosphere and waves are damped. However, the coefficient \(\kappa (\boldsymbol{x})\) will be set to zero in the area nearby the bodies so that damping will have no effect on the solution of the wave structure interaction problem. All the numerical models for the wave diffraction-radiation problem implemented in SeaFEM are based on [1] and [20]