Log In Sign Up

A spectral element solution of the 2D linearized potential flow radiation problem

by   Jens Visbech, et al.

We present a scalable 2D Galerkin spectral element method solution to the linearized potential flow radiation problem for wave induced forcing of a floating offshore structure. The pseudo-impulsive formulation of the problem is solved in the time-domain using a Gaussian displacement signal tailored to the discrete resolution. The added mass and damping coefficients are then obtained via Fourier transformation. The spectral element method is used to discretize the spatial fluid domain, whereas the classical explicit 4-stage 4th order Runge-Kutta scheme is employed for the temporal integration. Spectral convergence of the proposed model is established for both affine and curvilinear elements, and the computational effort is shown to scale with π’ͺ(N^p), with N begin the total number of grid points and p β‰ˆ 1. Temporal stability properties, caused by the spatial resolution, are considered to ensure a stable model. The solver is used to compute the hydrodynamic coefficients for several floating bodies and compare against known public benchmark results. The results are showing excellent agreement, ultimately validating the solver and emphasizing the geometrical flexibility and high accuracy and efficiency of the proposed solver strategy. Lastly, an extensive investigation of non-resolved energy from the pseudo-impulse is carried out to characterise the induced spurious oscillations of the free surface quantities leading to a verification of a proposal on how to efficiently and accurately calculate added mass and damping coefficients in pseudo-impulsive solvers.


page 1

page 2

page 3

page 4

βˆ™ 09/16/2020

Strong convergence of a Verlet integrator for the semi-linear stochastic wave equation

The full discretization of the semi-linear stochastic wave equation is c...
βˆ™ 06/22/2022

A monolithic Finite Element formulation for the hydroelastic analysis of Very Large Floating Structures

In this work we present a novel monolithic Finite Element Method (FEM) f...
βˆ™ 11/02/2019

A consistent discrete element method for quasi-static and dynamic elasto-plasticity

We propose a new discrete element method supporting general polyhedral m...
βˆ™ 09/11/2020

A time-spectral Stokes solver for simulation of time-periodic flows in complex geometries

Simulation of unsteady creeping flows in complex geometries has traditio...
βˆ™ 08/09/2021

Moment fitted cut spectral elements for explicit analysis of guided wave propagation

In this work, a method for the simulation of guided wave propagation in ...

1 Introduction

The offshore industry is of considerable economic and ecological importance in the world. In particular, floating offshore structures play a vital role in keeping up with society’s present and future demand for sustainable and green energy resources. Here will e.g. floating wind turbine foundations allow for the establishment of energy farms in greater water depths, where bottom-mounted foundations are too expensive to install [andersen2016floating]. Further development of wave energy converters, floating solar panel plants, and seaweed farms can also contribute in reducing the need for fossil-based fuels for energy production purposes [esmaeilzadeh2019shape].

1.1 Choosing the mathematical formulation

Designing structures at sea requires a very detailed, in-depth analysis of the physical environment to ensure a safe, well-functioning, and sustainable result. The study of waves and the interaction between waves and floating structures can be done in different ways, as it has been for many decades. One way is to use numerical methods and mathematical formulations that mimic the physical environment. The most precise mathematical formulation - yet computationally expensive - is the well-known Navier-Stokes Equations (NSE). The complexity and computational expense can be reduced dramatically by applying various assumptions to these equations, though, at a cost in terms of physical accuracy. A recent comprehensive review of wave energy converter design and different mathematical formulations is given in Davidson and Costello (2020) [davidson2020efficient].

The nonlinear and linear Potential Flow (PF) simplification is widely used for free surface gravity waves in the field of offshore and marine engineering. Compared to the NSE, the PF formulations will not be able to simulate breaking waves, overturning waves, nor viscous effects e.g. at the seabed or near structures; however, in offshore environments, it is often sufficiently accurate to model these effects through empirical boundary terms. The PF formulations require only a fraction of the computational effort needed by NSE-based solvers, and are thus much more suitable for long-time and large-scale wave simulations [ransley2019blind].

Considering wave-structure interaction and the linear PF setting, the classical impulsive time-domain formulation was widely used and studied in the 1980s and 90s [korsmeyer1988first, king1988seakeeping, bingham1994simulating]. King (1987) [king1987time] introduced the pseudo-impulsive formulation, and a significant amount of work on this topic has appeared over the past decade [read2011linear, read2012overset, read2012solving, amini2017solving, amini2018pseudo]. The basic concept of impulsive formulations is to force the floating body with an impulse in velocity, let the system evolve in time, and measure the resultant force on the body. A detailed investigation of the difference between the classical impulsive and pseudo-impulsive formulations is given in Section 5.2.

1.2 Choosing the numerical method for a potential flow formulation

Wave-structure interaction problems have traditionally been solved numerically using the Boundary Element Method (BEM) where the formulation is recast as an integral equation, which is solved - in a discrete sense - on the domain boundaries. A key feature of the BEM is the ability to easily resolve complex geometries. With the projection onto the boundaries, the formulation is reduced by one dimension. This essentially means that every boundary node receives information from the entire fluid domain, which leads to a dense non-symmetric system of equations. With this making large-scale simulation challenging due to the computational work effort when solved using conventional direct or iterative solvers. By using multipole expansions and/or pre-corrected FFT methods, the system can be made much more sparse, enabling asymptotic scaling, however with a large constant in front [jiang2012precorrected, white1994comparing]. The BEM has been implemented into the widely used state of the artfrequency-domain-based commercial software WAMIT [lee2006wamit]

; and also in the open-source software NEMOH

[babarit2015theoretical]. Impulsive time-domain implementations have as well been made using the BEM, e.g. for zero speed problems [korsmeyer1988first] and for forward speed problems [bingham1994simulating, korsmeyer1998forward, kara2011prediction].

Approaches for solving the PF problem have also been developed using the Finite Difference Method (FDM), both for pure nonlinear wave propagation [bingham2007accuracy, engsig2009efficient], nonlinear wave-structure interaction [ducrozet2010high], linear wave-structure interaction using pseudo-impulsive time domain formulations at zero speed [read2011linear, read2012solving, read2012overset], and with forward speed [amini2017solving, amini2018pseudo]. Using an FDM-based model, the entire fluid domain is discretized to solve the Laplace problem, hereby resulting in a larger linear system of equations compared to the BEM discretization; however, the system is sparse due to the locally applied stencils. This can - via iterative solvers - lead to an optimal scaling in work effort for higher order accurate schemes [bingham2007accuracy, engsig2009efficient]. Complex geometries can also be represented using this discretization method; yet their inclusion is far from trivial since updating and adapting the mesh domain and imposing boundary conditions consistently can be challenging within this discretization framework.

The Finite Element Method (FEM) has also been used for modeling linear and nonlinear wave propagation and wave-structure interaction [wu1994finite, taylor1996analysis, ma2001finite, wu2003coupled]. Due to the local support of the linear basis functions used in the FEM, the resulting linear system of equations will be sparse, hereby enabling optimal computational scaling and a small memory requirement. Also, the element discretization approach makes it very suitable for complex geometries; however, re-meshing strategies may be necessary in the nonlinear setting due to the movement of the free surface boundary. A drawback for the FEM is the restriction to second order in spatial accuracy from the use of linear basis functions.

This caveat naturally leads to another element based method, namely the higher order version - or extension - of the FEM, i.e. the spectral element method (SEM). This method can also be seen as a multi-domain version of the classical single domain spectral method [canuto2007spectral, canuto2007spectral_evolution, kopriva2009implementing], where higher order polynomial basis functions are used to enable higher order spatial accuracy. The use of SEM for fluid dynamics problems was originally proposed by Patera (1984) [patera1984spectral]. In the late 90s, the SEM was shown to be inherently unstable for free surface flows using simplex elements by Robertson and Sherwin (1999) [robertson1999free]; however, almost two decades later, the first stable SEM-based fully nonlinear PF model for pure wave propagation was successfully developed by Engsig-Karup et al. (2016) [engsig2016stabilised]. Since this breakthrough, multiple extensions have been developed and applied to model wave-wave, wave-bed, and wave-structure interaction with optimal scaling through -multigrid methods [engsig2019mixed, laskowski2019modelling, engsig2021efficient, mortensen2021simulation]. A recent review of the state of the art for SEM models can be found in Xu et al. (2018) [xu2018spectral]. The clear advantage of the SEM is the geometrical flexibility, optimal computational scaling properties, and high spatial accuracy. In combination with the general increase in computational resources, long-time and large-scale simulations can be achieved using the SEM as the numerical discretization method for linear and nonlinear PF formulations.

1.3 Paper contributions and overview

In this work, we develop a linear potential flow spectral element (PF-SEM) solver for the simulation of linear wave-structure interaction through a pseudo-impulsive PF formulation. The unstructured mesh model has the ability to take into account the geometric complexity of the structures by using affine and curvilinear elements and is validated through convergence studies, stability considerations, and simulations that are compared against known benchmark results from the literature and state of the art solvers. The computational effort is shown to scale with , where . This confirms the potential of the model and completes and extends the preliminary results presented by Visbech et al. (2022) [visbech2022efficient]. An additional novelty is the investigation of the non-physical behavior of the numerical solution at the free-surface/body intersection arising from unresolvable energy in the pseudo-impulse (or impulse), highlighting the importance of properly tuning the impulse signal.

The rest of the paper is organised as follows: In Section 2, we present the mathematical formulation in terms of the governing equations for the pseudo-impulsive formulation. Then, in Section 3, we outline the numerical method and discretization, with this covering the SEM, temporal integration, computation of pressure and hydrodynamic coefficients, and the construction of the discrete pseudo-impulse. In Section 4 the numerical properties are outlined by showing studies of convergence and computational scaling effort, but also, highlighting linear stability analysis of the spatial discretization. Section 5 seeks to display different numerical experiments obtained using the proposed model, and in Section 6, a conclusion is made to summarize the novelty achievements and findings.

2 Mathematical Formulation

We consider a spatial 2D fluid domain, , upon which the flow is assumed to be: 1) impressible, 2) inviscid, and 3) irrotational. These simplifications enable the NSE to be represented via PF with the fluid velocity field, , expressed as the gradient of a scalar velocity potential, , where is the spatial differential operator

in a Cartesian coordinate system spanned by the horizontal

-axis and the vertical -axis. Let denote the time domain by , with being the time. Now and .

An even further simplification of the PF equations is often used in offshore engineering by assuming that the waves and the oscillatory motions of the floating structure are of small amplitude compared to the wavelength and the size of the structure, respectively. This yields the linearized version of the PF equations, where the perturbation expansion term of first order of the complete nonlinear solution still captures the majority of the physics. The linear assumption enables the general solution for wave-structure interaction to be decomposed into independently solvable unsteady radiation and diffraction parts [haskind1946hydrodynamic], where the focus of this paper is on the former. Thus


with denoting the incident wave potential and the scattered wave potential caused by the interaction of the fixed structure with the incident waves. Finally, denotes the total radiation potential,

, due to rigid-body motions in six degrees of freedom, i.e. surge (

), sway (), heave (), roll (), pitch (), and yaw (). For the 2D setting only surge, heave, and pitch are relevant.

Figure 1: Illustration of the 2D fluid domain, , and the associated boundaries, . The free surface is denoted by and the water depth by . The arbitrary domain length in is denoted by .

The entire fluid domain, , will be bounded by five different boundaries, , as illustrated on Figure 1, namely: 1) from above at is the time-dependent free surface boundary, , 2) from below at by the bathymetry boundary, , 3) from the ends of the domain by the far-field boundaries, , 4) from the floating body by the body boundary, , and in some special instances 5) from the symmetry boundary, . We have that is the free surface elevation and is the water depth. For the entirety of this work, the depth is assumed constant, .

2.1 Governing equations

The governing equations for each of the linear radiation potentials, , is achieved using the pseudo-impulsive formulation by combining the linear Eulerian kinematic and dynamic free surface boundary conditions


with the mass conserving Laplace boundary value problem including suitable boundary conditions as


where m/s is the constant gravitational acceleration, is the derivative in the normal direction to the boundary, is the pseudo-impulsive velocity from the Gaussian displacement, , that will be defined in Section 3.4, and

is a generalized unit normal vector given as


where is the unit normal vector and is the position vector of a point on .

To complete the Laplace problem and make it well-posed, boundary conditions on the symmetry boundary, , and the far-field boundary, , are needed. If is applicable, it depends on the solution of the problem as


where the former applies for and , and the latter for . On the main object is to ensure that no energy is reflected back into the computational domain, ultimately disturbing the solution. For this work, two different approaches are used to make a non-reflective domain: 1) An impermeability condition, , on combined with a horizontal grid stretching strategy [xxx], or 2) a classical relaxation zone [larsen1983open] (for the high frequency content) combined with a Sommerfeld condition [zhang2009non], , on (for the low frequency content). Assuming a straight and vertical far-field boundary, the imposed velocity, , at time is determined as , with being determined from the shallow water limit. The discrete time step size, , is defined in Section 3.2.

2.2 The pseudo-impulsive radiation problem

From each of the computed radiation velocity potentials, , the dynamic pressure, , where , can be found from the linearized Bernoulli equation at every point in for any time instance in . From this, the radiation forces, , can be evaluated by integrating over the wetted-body surface, . With this, the radiation force in the ’th direction due to a pseudo-impulsive motion in the ’th direction is


where is the generalized normal vector given by (7).

Let denote the frequency space , with being the radial frequency that is connected to the cyclic frequency as . Following Amini-Afshar and Bingham (2017) [amini2017solving], the frequency-dependent added mass and damping coefficients, , are given by


where denotes the Fourier transform operator, .

3 Numerical Methods and Discretization

Exploiting a traditional Method Of Lines (MOL) approach, the governing equations are discretized spatially and kept continuous in time. We seek to represent the spatial part by the SEM, ultimately yielding a semi-discrete system of ordinary differential equations, which then are time-integrated using the classical Explicit 4-stage 4th order Runge-Kutta Method (ERK4). The time-derivative for the pressure computation needed in (

9) will be performed by a 4th order finite difference approximation, and the Fourier transformations in (10) are handled discretely using the Fast Fourier Transform (FFT)[cooley1965algorithm]. Finally, the design of the discrete pseudo-impulsive velocity and displacement used in (6) and (10), respectively, will be specified. All of the aforementioned objectives are outlined in the following.

3.1 Spatial discretization - The spectral element method

The fluid domain, , will be tessellated using conforming non-overlapping affine or curvilinear triangular elements in an unstructured manner, where we can arrange the elements uniformly or apply mesh refinement at the free surface and body boundary, and . See Figures 2, 3, and 7 for visualizations. The tessellation will consist of elements such that we, at some particular time instance, have a time-invariant partitioning of the domain and the global radiation velocity potential as


with being the ’th element having sets of coordinates , and is the local solution defined on this element. We seek to apply the following approximate polynomial representation of the local solution


where the former and latter sum expressions are the modal and nodal forms, respectively. is the number of local expansion coefficients, , or local solution values, . Moreover, is the modal basis functions (Jacobi polynomials) up to order and is the nodal basis functions (Lagrange polynomials) up to order . The order of the basis functions, , is connected to as . Each element is mapped onto a 2D simplex reference domain, , given by: . Affine or curvilinear mappings, and , are introduced, ultimately giving rise to a local transformation Jacobian, . On , an orthonormal basis, , is constructed on a nodal set of grid points, , as outlined in Hesthaven and Warburton (2007) [hesthaven2007nodal]

, using the family of Jacobi polynomials and Gauss-Lobatto points. Hereby leading to a well-behaved generalized Vandermonde matrix for interpolation


3.1.1 Weak formulation and Galerkin spectral element discretization

An essential step in the quest for the spatial discretization of the pseudo-impulsive formulation using the SEM is the derivation of the weak formulation in integral form. Ultimately, this results in the removal of second-order derivatives and enables the natural incorporation of Neumann boundary conditions into the scheme. Dividing the boundaries into Dirichlet and Neumann boundaries as , and considering a test function on as , that is assumed to be infinitely smooth and vanishing on . By: i) multiplying the time-invariant Laplace problem in (3) by , ii) applying integration by parts, iii) using the properties of on , and iv) rearranging, the weak formulation yields


where is the Neumann flux which is only relevant for inhomogeneous conditions on and in the case of a Sommerfeld condition. A Galerkin approach is adopted by setting the test function equal to the basis function, , in the weak formulation, and then inserting the polynomial representation from (12). With this, a global linear system of equations of Degrees Of Freedom (DOF) is obtained in the form


here is the sparse system matrix, is the system vector containing information about the boundary conditions, and being the number of nodal points in the spatial discretization. Local contributions to the system matrix, , can be derived[engsig2016stabilised] as


with being the local mass integration matrix, , and and is the differentiation matrices in the - and

-direction, respectively, defined through the chain rule as

and . Now, the differentiation operators on are given as and with and for . Similar contributions - this time in a 1D setting due to the line integral - can be derived for the system vector as


where is initialized to be all zeros, is the ’th mass integration matrix on the boundary element, and are the discrete nodal flux values.

Inhomogeneous Dirichlet boundary conditions on are not enforced due to the vanishing properties of the test function, , on . These are strictly imposed by modification of the system matrix and vector, and , in standard ways.

3.1.2 Curvilinear elements and non-constant transformation Jacobians

For curved physical boundaries a geometrical second order error is introduced when using affine elements; ultimately, ruining the favorable spatial spectral convergence rate of the SEM. This challenge can be solved by considering curvilinear elements that will be constructed by modification of existing affine elements through transfinite interpolation with linear blending due to Gordon and Hall (1973) [gordon1973transfinite, gordon1973construction].

The curvilinear elements will - in contrast to the affine elements - result in non-constant transformation Jacobians leading to the possibility of aliasing errors in the approximation due to the nonlinear transformation in the discrete inner products, and as a result of this affecting the convergence rate. To overcome this challenge, a super-collocation method [kirby2003aliasing, engsig2016stabilised, engsig2019mixed] is incorporated to handle the discrete inner products.

3.2 Temporal integration

The method of lines approach leaves us with a semi-discrete system of ordinary differential equations from (2), that are evolved in time using the classical ERK4. In combination with the time integration scheme, a global Courant Friedrichs Levy (CFL) condition is exploited to ensure sufficient stability, accuracy, and efficiency. Hereby, an upper bound for the time step size, , is set to be


where is the Courant number, is the minimum grid spacing on the free surface, and is set to be the asymptotic shallow water limit. This choice of time step size can be seen to be quite conservative, as the propagation of higher frequency waves will need multiple time steps to move an equivalent . Using this time step size, the time domain, , is partitioned uniformly by sample points


3.3 Computation of the pressure and hydrodynamic coefficients

As is evident from (9), the computation of the radiation force, , require information about the dynamic pressure, . The evaluation of this quantity involves time derivatives of the radiation potentials that will be computed discretely by finite difference approximations [leveque2007finite] as


where the time vector, , is discretized uniformly as in (19) given the time step size from (18). The coefficients are determined though the method of undetermined coefficients dependent of the order of accuracy and whether forward, centered, or backward stencils are considered. To complement the th order accurate time integration scheme, an equivalent accurate finite difference scheme is constructed using different stencils, and the convergence of the derivative approximation has been successfully validated.

The computation of the Fourier transforms needed in (10) for the evaluation of the hydrodynamic coefficients is applied discretely through FFTs [cooley1965algorithm]. The frequency domain, , is discretized uniformly with the frequency resolution, , with being the time length of the simulation. To increase the resolution, the displacement and force signals, and , are extended with zeros assuming that the latter have reached a zero steady state solution, which will be ensured by following the design proposal in the upcoming section. This particular technique has been used extensively throughout the presented work.

3.4 Construction of the discrete pseudo-impulsive displacement

On the body boundary, , a pseudo-impulsive velocity, , is to be imposed. This velocity will be based on a Gaussian displacement signal, , [amini2017solving, amini2018pseudo] and designed such that: 1) the Gaussian is tailored to contain a specific range of frequencies, 2) the signal is practically zero at , and 3) the resultant radiation force signal, , should reach a zero steady state solution at . With these requirements, a unit height Gaussian displacement signal is defined as


where is the time location of the Gaussian peak and is a frequency range governing parameter. The parameter, , is some small quantity to ensure that and is governing how many orders of magnitude there should be between the frequency components and , where is the chosen maximum resolvable design frequency so that the Gaussian contains frequencies between and . For this work, has been used. The maximum resolvable design resolution is defined through linear wave theory as


with being the largest distance between two nodal points on the free surface boundary, , and is a parameter governing the length of the smallest resolvable design wavelength. Unless stated otherwise, , resulting in grid points per smallest possible wavelength. Analytical expressions for the pseudo-impulsive velocity and the Fourier transform of the displacement signal can be obtained as:


Finally, the simulation time is set to be at least to ensure that .

4 Numerical Properties

In the following, the numerical properties of the developed model are highlighted. This will be done by considering: 1) the spatial accuracy of the Laplace solver in terms of a convergence study, 2) the efficiency when solving the sparse linear system of equations in (15), and 3) the temporal stability in relation to the spatial discretization.

4.1 Convergence of Laplace solver

As an exact solution to the pseudo-impulsive radiation problem is non-trivial, we consider the convergence study of the Laplace solver in the scope of the Method of Manufactured Solutions (MMS) [roache2002code], where a true solution, , is assumed to the problem, from which analytical expressions can be derived for the boundary conditions and right hand side function, and hereby reforming the Laplace problem into a Poisson problem. The true solution should be infinitely smooth, e.g. a combination of trigonometric functions, to complement the wave problem, and with this, we seek to verify that we are able to archive: 1) spectral -convergence for fixed meshes by increasing the order of orthonormal basis functions, , and 2) algebraic -convergence of order for fixed orders of orthonormal basis functions with decreasing element size.

For the convergence study, we consider a fluid domain, , that - due to the symmetry boundary - has a quarter of a floating cylinder with radius, , in the top left corner. The domain is partitioned using curvilinear elements around to capture the curved boundary of the structure, and for the remaining part, we seek to use regular affine elements in an unstructured fashion. When considering the convergence study, the elements are arranged uniformly, as seen in Figure 2 for the three coarsest meshes used having a total of number of elements, respectively. In Figure 3, a zoomed-in plot of the body boundary elements for the mesh with is pictured to show how the distribution of nodal points are placed within the curvilinear elements for .

Figure 2: Three unstructured uniform meshes using triangular elements with decreasing element size, , for discretization of a floating cylinder with radius .
Figure 3: Nodal distribution (red dots) of affine and curvilinear elements near for different polynomial orders, , on the coarsest mesh in Figure 2, .

The convergence study is performed as shown in Figure 4 and Figure 5 with the - and

-test, respectively, where the error is evaluated globally in the solution to the partial differential equation using the infinity norm as

. For the former test, the three meshes in Figure 2 are used with increasing order of the orthonormal basis functions, . The spectral convergence behavior is confirmed in the semi-logarithmic plot by comparing with a line, where the error can be seen to reach about , from where round-off errors start to become significant and accumulate. It can also be noted that the finer meshes, in general, give smaller errors than the coarser ones. This is an early indication of -convergence, which ultimately can be confirmed in the latter test on the logarithmic plot, showing algebraic convergence of order for by decreasing the size of the elements.

Figure 4: -convergence study using curvilinear and affine unstructured triangular elements. Fixed meshes from Figure 2 with increasing polynomial order, .
Figure 5: -convergence study using curvilinear and affine unstructured triangular elements. Fixed polynomial order, , and decreasing elements size in terms of maximal element edge length, .

4.2 Scaling of the computational effort

When solving free surface potential flow, the main concern and computational bottleneck is the work associated with solving the Laplace problem numerically, which has to be done multiple times - dependent on the time discretization method - at every time step. Therefore, one should seek to complete this challenge as efficiently as possible.

The linear system of equation in (15) is solved directly in MATLAB, , using a sparse reverse Cuthill-McKee reordering [CuthillMckee]. The permutation, , reorders the system matrix as , such that the band width is minimized, hence making the direct solver much more efficient. In Figure 6, the study of computational effort can be seen in terms of the DOF of the system, , against CPU time per time step. Nine different uniform meshes where constructed with the ratio between the coarsest and the finest mesh as yielding meshes having number of elements. For each mesh, polynomial orders of were applied, resulting in . From the study, can be observed with , when solving the linear system using a direct solver. The reason why is due to fill-in of by the direct solver.

Despite the fact that the sparse scheme - due to the local support of the basis functions - already is very efficient, the computational effort can be made to scale with by considering the use of a geometric -multigrid method that exploits the polynomial convergence features of the SEM model to solve the Laplace problem iteratively [laskowski2019modelling, engsig2021efficient].

Figure 6: Analysis of computational work effort for the 2D SEM solver using a sparse reverse Cuthill-McKee ordering. Each simulation (black dot) represent an average of solving the Laplace problem times.

4.3 Linear stability analysis of spatial discretization

Temporal instability is not solely related to temporal discretization as outlined in Robertson and Sherwin (1999) [robertson1999free] where the issues related to computational instabilities of incompressible and inviscid free surface flow are discussed in the context of the SEM. The reason was found to be caused by asymmetrical spatial discretizations of the free surface of the fluid domain, for which the authors proposed a method for studying the stability of the scheme. Due to the adopted MOL approach, we are left with a semi-discrete system for the kinematic and dynamic free surface condition from (2

). The eigenvalues and eigenvectors of this system will dictate the temporal stability of the scheme, where eigenvalues with positive real components indicate an unstable system, as these govern the diffusive properties, thus causing the system to diverge. Such analysis is performed in a recent work


One way to ensure stability when solving free surface flows using the SEM is by applying quadrilateral elements for the entirety of the domain [engsig2016stabilised, mortensen2021simulation]; however, in Engsig-Karup et al. (2019) [engsig2019mixed] it was shown to be sufficient to exploit a single layer of this type of elements below the free surface boundary, . It was found that the instability issue is related to the accuracy of computing the vertical free surface velocity, which govern the dispersive property of the discrete model. To maximize this, the nodal distribution beneath must be on a vertical line. In the same paper, weak stability properties, when using triangular elements, were highlighted, showing that some polynomial orders were keen to provide instability issues despite the fulfilled symmetry requirement. Similar analyses have been carried out in relation to this work, and similar results have been obtained. Hereby emphasizing the fact that stability analyses are solely related to the specific meshes; with this, one should be careful - when using triangular elements - picking the proper polynomial order that allows for a stable temporal evolution of the governing free surface equations.

5 Numerical Experiments

We seek to outline simulation results for floating bodies (a cylinder and a box) in the following by considering comparisons of force signals and added mass and damping coefficients with other well-established numerical models and benchmark results from the literature. Also, an extensive investigation of the numerical phenomenon of free surface spurious oscillations arising from unresolved energy in the pseudo-impulse is carried out.

5.1 Heaving and surging structures

For the simulations presented in the following, affine and curvilinear elements have been used, where the elements on the free surface boundary, , is constructed symmetrically to stabilize most of the polynomial orders. Only stable orders have been considered. The elements will be refined toward the free surface and body boundaries, and , e.g. as shown in Figure 7. On the far-field boundary, , a combination of a Sommerfeld condition [zhang2009non] and an absorption zone technique [larsen1983open] are used to make the domain non-reflective for this multi-frequency wave problem.

Figure 7: Visualization of mesh used for simulation of a floating cylinder. Mesh refinement on and with 5 curvilinear elements constituting for meshing of a quarter of the cylinder due to symmetry in the problem. The radius is m, the domain length is m, and the finite water depth is m.

5.1.1 The floating cylinder

The first floating structure is the classical half submerged cylinder with a radius . We compare the calculations with an equivalent FDM-based solver, used to obtain the results in Read and Bingham (2012) [read2012overset] and others [read2011linear, read2012solving], by exposing the two models to identical input displacement signals, , and then studying the corresponding output force signals, . For the simulation, we used the mesh shown in Figure 7, having 5 curvilinear elements around for meshing of a quarter of the cylinder surface. Basis functions of order are chosen to ensure stability. From the comparison in Figure 8, an excellent visual agreement can be observed as the two force signals are placed on top of each other yielding non-dimensional absolute errors in the order of . The minor discrepancies are located around the peaks of the force signals and are associated with differences in spatial resolution.

Figure 8: Comparisons of heave-heave radiation force signals, , of the FDM-based and the SEM models. Simulation obtained using ordered basis functions on the mesh shown in Figure 7. The design wave period for the displacement signal is denoted by .

Next, a comparison of added mass and damping coefficients with results published in Read and Bingham (2012) [read2012overset]. These benchmark results have been validated against analytical infinite-depth solution calculated using a multipole method [ursell1949heaving]. The radius of the cylinder is set to m and the water depth is m yielding [-]. The added mass and damping coefficients have been normalized as


The surge-surge and heave-heave results are displayed in Figure 9 (left) and 9 (right), respectively. The numerical results show close to perfect visual agreement, especially for where the waves - in terms of wave length - are sufficiently short to be uninfluenced by the finite depth. The latter comparison indeed highlights the legitimacy of the proposed model.

(a) Surge-Surge case: and .
(b) Heave-Heave case: and .
Figure 9: Dimensionless added mass and damping coefficients for a floating cylinder of radius m and a finite water depth of m. Basis functions of order is used. Comparison between FDM-based solver results [read2012overset, read2012solving] and SEM model.

5.1.2 The floating box

The second floating structure to consider is a floating box-type structure with the draft, , and the half-length, . A comparison with the analytical method in an infinite fluid domain of finite water depth presented by Zheng et al. (2004) [zheng2004radiation] is to be carried out. The aforementioned method was validated against Lee (2005) [lee1995heave] and a BEM model.

For the simulation, we adopt the same setting as in the comparison study, hereby using a dimensionless depth of [-] and dimensionless length of [-]. As all boundaries are non-curved, solely affine elements will be considered. Regarding the meshing of the box, 6 elements constitutes the bottom of the structure and 12 elements the side. Basis functions of order is used. The added mass and damping coefficients have been normalized as


In Figure 10 (left) and (right) the numerical results for non-dimensional added mass and damping coefficients can be seen for surge-surge and heave-heave, respectively. For , an excellent visual agreement can be discovered as for the cylinder case. Note that the results from Zheng et al. (2004) have been manually digitized.

(a) Surge-Surge case: and .
(b) Heave-Heave case: and .
Figure 10: Dimensionless added mass and damping coefficients for a floating box with dimensionless depth, [-], and dimensionless length, [-]. Basis functions of order is used. Comparison between Zheng et al. (2004) [zheng2004radiation] and the SEM model.

5.2 Investigation of spurious oscillations

As presented in the introduction, the difference between the classical impulsive formulation [korsmeyer1988first, king1988seakeeping, bingham1994simulating] and the pseudo-impulsive formulation [read2011linear, read2012overset, read2012solving, amini2017solving, amini2018pseudo] will be highlighted in the following. Starting with the former, we have that the impulsive motion is imposed rapidly at initial time, , using the generalized Dirac Delta function, , i.e. a pressure release. From the Fourier transformation, we conclude that , hence unit amplitude energy is imposed on all frequencies in the spectrum. This operation works well in the continuous setting; however, when the governing formulation is discretized, one faces the challenge of resolving energy on frequency components that are unresolvable by the discrete scheme. From experiments with the classical impulsive formulation, the unresolved energy has been seen be alias towards lower frequencies and give rise to spurious oscillations of the free surface quantities.

This is the primary motivation for using the pseudo-impulsive formulation, where the frequency spectrum is tailored through a unit-amplitude Gaussian displacement function as outlined in Section 3.4. From the experiences with the classical formulation, we seek to study the unresolved energy in terms of the spurious free surface oscillations using the pseudo-impulsive formulation as this - to the authors’ knowledge - has not been investigated before. The study will be carried out by including too high frequencies; that will enable us to explore the behavior of the system. For the analysis, we solely consider the case of a floating cylinder in a heave-heave motion and compare it with a BEM-based model on infinite depth. The BEM-model has been validated against solutions on infinite-depth computed using the method of Porter (2008) [porter2008solution]. Also, a stable polynomial order of is used on a mesh with - until further notice - curvilinear elements around a quarter of the cylinder surface. The finite water depth is set to [-].

5.2.1 Pseudo-impulse response

We first investigate the response when the same discrete system is exposed to different pseudo-impulses. Recall the parameter, , from (22) governing the smallest resolvable wavelength, , where is the largest distance between two grid points on the free surface. For the results presented so far, a value of has been used with the argument that a linear wave requires (approximately ) points to be represented adequately. Thus, we arrive at . In the following, we change the -parameter in the range , yielding both under- and over-resolved energy spectra. All of the displacement signals, , can be seen on Figure 11 (left), where the same peak location, s, have been used. The corresponding force signals, , can be seen on Figure 11 (right). Observe how decreasing implies larger gradients in the displacement signal, hence larger boundary imposed velocities, ultimately yielding higher forces. Recall that a narrow displacement signal implies a wide band of frequencies, and vice versa for a wide signal.

Figure 11: Left: Dimensionless displacement signals, , for . Right: Dimensionless force signals, , for .
Figure 12: Top left: Dimensionless free surface elevation, , at the point between and for . Top right: Zoomed-in snapshot of the top left plot. Bottom: Fourier transformation of for .

When studying the free surface elevation, , right at the point between and , , the consequence of using different pseudo-impulses can be seen in Figure 12 (top left and right). Here the presence of the spurious oscillations starts to take place for , which confirms the initial use of to match the smallest resolvable wavelength. Also, the occurrence of the oscillations can be seen to be of constant phase independent of the pseudo-impulse; yet, the amplitudes tend to increase for decreasing . On Figure 12 (bottom), the Fourier transformation of the free surface signal, , is shown, ultimately confirming the aforementioned observation of same phase and amplitude when decreases as seen in increased spectral energy at the same particular frequency, . With this, it can be concluded that the spatial discretization governs which frequency the unresolved velocity imposed energy is aliased towards. Furthermore, it can again be verified that seems to provide the right physical results. Remark: Identical behavior have been observed in terms of the free surface velocity potential, .

On Figure 13, the dimensionless added mass and damping coefficients, and , are plotted for . It can be observed how narrowing the displacement signal yields an increased frequency span. Note that the results for goes out to , hence not shown on the figure. In general, it is evident that near the cut-off frequency, , some other oscillations are present. These are caused by the fact that the added mass and damping coefficient are computed via: . Hence near , a small number is divided by another small number, thus prone to give rise to this kind of disturbance. Another result is how the added mass and damping coefficients can be well-approximated by decreasing the value of - even to a level with the smallest introduced waves being far from adequately resolved by the spatial scheme. Somehow, the non-physical numerical time-domain spurious oscillations do not impact the frequency solution. The apparent advantage of this observation is that one can evaluate added mass and damping coefficients extremely efficiently by choosing sufficiently small values of , resulting in a narrow displacement signal, requiring less simulation time and a smaller domain size. Remark: The efficient computational approach was tried for the FDM-based pseudo-impulsive model [read2011linear, read2012overset, read2012solving], and similar result was discovered.

Figure 13: Top: Dimensionless added mass coefficients, , for . Bottom: Dimensionless damping coefficients, , for . Both compared with an infinite-depth BEM model.

5.2.2 Spatial resolution

The next to investigate is the dependency on spatial resolution. Choosing the displacement signal to be very narrow by setting in (21) will enforce oscillations of the free surface quantities. Then, we seek to apply different meshes with increasingly more elements at and around . Let the parameter, , govern the number of elements on a quarter of the cylinder, where the element edge length on is (approximately) the same as on . Recall from the previous investigation that was used, whereas is considered in the following.

Due to the identical displacement signals, the force signals are - to a large extend - identical; yet, minor discrepancies will occur due to resolution. When considering the free surface elevation, , shown in Figure 14 (top left and right), the spurious oscillations are very much present for all simulations; however, the phases are not equal for different , and the amplitude can be seen to decrease for increasing . Also, the frequency seems to increase for increasing . In Figure 14 (bottom) the equivalent frequency content of the free surface elevation, , can be seen. Here, the time-domain observations can be confirmed, as it is clear that by increasing the spatial resolution: 1) the energy of the spurious oscillation is decreased, and 2) the location of the oscillation frequency component migrates towards higher frequencies. With this, it can be concluded that by refining the spatial resolution in the presence of unresolved energy, the frequency, on which the energy is aliasing towards, becomes higher and less significant.

Figure 14: Top left: Visualization of free surface elevation, , at the point between and for . Top right: Zoomed-in snapshot of the top left plot. Bottom: Fourier transformation of for .

6 Conclusion

A novel 2D spectral element model is proposed for the simulation of linear radiation potentials in relation to floating offshore structures. The mathematical formulation was established within the framework of linear potential flow using a pseudo-impulsive time domain formulation. The numerical methods and solution approaches were established for solving the governing equations numerically, where the spectral element method is handling the discretization of the spatial domain.

The numerical properties of the model were tested in terms of spatial accuracy, computational efficiency, and temporal stability, with this confirming - and -convergence, scalability, where , and instability considerations, respectively. At last, various numerical experiments were carried out and compared against established solvers and benchmark results from the literature. All showing excellent visual agreements and ultimately confirming the legitimacy of the proposed model. Also, a thorough novelty investigation of non-resolved velocity imposed energy leading to spurious oscillations of the free surface quantities was performed. With this, concluding and emphasizing the importance of correctly tailoring the Gaussian input signal, and also discovering an efficient non-physical way of computing added mass and damping coefficients.