A high-order spectral element unified boussinesq model for floating point absorbers

Nonlinear wave-body problems are important in renewable energy, especially in case of wave energy converters op- erating in the near-shore region. In this paper we simulate nonlinear interaction between waves and truncated bodies using an e ﬃ cient spectral / hp element depth-integrated uniﬁed Boussinesq model. The uniﬁed Boussinesq model treats also the ﬂuid below the body in a depth-integrated approach. We illustrate the versatility of the model by predicting the reﬂection and transmission of solitary waves passing truncated bodies. We also use the model to simulate the motion of a latched heaving box. In both cases the uniﬁed Boussinesq model show acceptable agreement with CFD results – if applied within the underlying assumptions of dispersion and nonlinearity – but with a signiﬁcant reduction in computational e ﬀ ort.


INTRODUCTION
This work is concerned with the development of an efficient tool for nonlinear analysis of floating wave energy converters (WECs). Floating WECs are structures that harvest energy from the wave motion. There is a multitude of different designs and operation modes of WECs, see e.g. Greaves and Iglesias (2018) for an overview. In particular, we are interested in developing methods for point absorbers, i.e. axi-symmetric devices having a small diameter compared to the characteristic wave length in which they operate.
The tools-of-the-trade for predicting the motion and power production of WECs are models solving the Cummins equation (Cummins, 1962) using hydrodynamic coefficients computed from linear potential flow. These models, however, are fundamentally based on the small-amplitude/small-motion assumptions. Thus the models struggle with nonlinear effects important for survival cases, varying bathemetry and WECs operating under phase control (or simply operating in the natural resonance region). During the last, say, five years, Reynolds-Averaged Navier-Stokes (RANS) simulations have become popular for WECs, e.g. (Yu and Li, 2013;Palm et al., 2016;Ransley et al., 2017). RANS is a complete and accurate model with respect to nonlinear phenomena but computationally very costly. For example, simulations of 3 hour irregular sea states have been reported to take in the range of 30 000 − 150 000 CPU hours (Eskilsson et al., 2015b;Kim et al., 2016).
In this paper we outline a body-wave interaction model based on depth-integrated Boussinesq-type wave equations (Peregrine, 1967;Abbott et al., 1978;Madsen and Sørensen, 1992;Madsen and Schäffer, 1999;Brocchini, 2013). Numerical models based on Boussinesq-type equations are standard engineering tools for simulating nonlinear wave propagation and transformation in coastal areas. They are popular because they are computationally efficient as the equations are expressed in horizontal dimensions only. The depth-integration also circumvent the problem of a time-dependent computational domain caused by the moving free surface boundary condition. To include truncated bodies in the the depth-integrated approach is not straightforward. Typically floating bodies are introduced in some simplified manner (Ertekin et al., 1986;Chen and Sharma, 1995;Nørgaard and Andersen, 2012).
Our study is based on the work of Jiang (2001) on the 'unified' Boussinesq model. Jiang decomposed the domain into a free-surface domain and a body domain. The term 'unified' comes from the fact that Jiang modelled also the domain under the body with a depth-integrated approach. A similar setting was considered by Lannes (2017a). Lannes study mainly dealt with the hydrostatic nonlinear shallow water equations and deduced semi-analytic nonlinear solutions for the wave-body problem.
In this paper we outline the governing equations and the resulting numerical scheme. We then apply the model to the case of a solitary wave impacting a pontoon (Lin, 2006) and investigate the effect of different pontoon shapes. Additionally, we use the model to simulate phase control of a heaving WEC. Latching (Budal et al., 1979) makes the velocity of the body to be in phase with the wave excitation force, increasing the response and possible power production.

HYDRODYNAMIC MODEL
Shallow water and Boussinesq models for free surface flows can be derived from the fully nonlinear potential equations for an incompressible, irrotational and non-viscous fluid by expanding the velocity potential in terms of the vertical coordinate and integrating the Laplace equation over the water depth. The result is a set of equations expressed in horizontal dimensions only. Please note that there exists a plethora of depth-integrated models with different assumptions on the nonlinearity and dispersive properties, see e.g. Madsen and Schäffer (1999).
The unified approach is based on domain decomposition. We divide the total domain into an outer free surface sub-domain (Ω w ) and an inner sub-domain (Ω b ) that represents the area under the structure, as shown in figure 1. The model equations, valid in both domains, for the case of constant bathymetry read (Jiang, 2001;Lannes, 2017a) where the variables d(x, t), u(x, t) and q(x, t) are the depth measured as the height of the water column, depth-averaged velocity, and the mass flux (q ≡ du), respectively. We have denoted with g the acceleration of gravity, and by P = gd + Π the total pressure. Here Π is the pressure at the free surface in Ω w and on the body bottom in Ω b . In the free surface region Π is equal to the atmospheric pressure and set to zero, but on the body surface the pressure is neither constant nor known. Hence The term Φ take into account the dispersive terms in the two domains. Please note that for a purely heaving body there is no higher order terms present so the equation simplifies to the hydro-static nonlinear shallow water equations (NSW) in the body domain (Lannes, 2017a). In the free surface domain we use the Madsen and Sørensen (1992) (MS) formulation. Subsequently, The tunable parameters are defined in the literature as α MS = 1/15 and B = 1/3 + α MS and will optimize the linear dispersion relation of the system (Madsen and Sørensen, 1992). The pressure P(x, t) in the inner/body domain can be evaluated by taking the divergence of eq. (1b), and combine with the time derivative of eq. (1a). Assuming only heave motion we introduce the vertical acceleration a = d tt . Thus, we can show that the total pressure satisfies the following equation The main difference between the two domains is that in Ω w the total pressure and the free surface elevation are readily obtained by eqs. (1a)-(1b), automatically satisfying eq. (4). On the other hand, in the inner domain Ω b , the relation (1a) acts as a constraint on the flux divergence, exactly as in incompressible flow. The total pressure is in this case determined from eq.s (1b) and (4). For a purely heaving body, the vertical acceleration will be determined by the application of Newton's second law to the body where m b the mass of the body and ρ w is the water density. Unfortunately, the transmission conditions between the dispersive outer domain and the non-dispersive inner domain cannot be rigorously formulated in the nonlinear case (Lannes, 2017a,b) and we have decided to handle them numerically. We do this by introducing a small transition layer Ω l on the sides of the body, where we solve the NSW equation. The size of the layer is calibrated to be long enough such that we avoid the propagation of dispersive terms under the body, where they are equal to zero, and short enough to permit the propagation of the wave with minimal or no distortion.

NUMERICAL MODEL
The focus of our study is in 2D (vertical plane) waves modelled by a 1D system of PDEs. In the outer domain, we will solve the 1D MS model expressed in a first order formulation where we have multiplied the mass equation (6) by g such that we can use the same set of variables (P, q), through all the domains. In the inner domain we solve the non dispersive 1D NSW system In each domain, we employ the continuous spectral/hp element method (Karniadakis and Sherwin, 2013) to solve the systems (6) and (7). We describe the discretization process only to one equation, knowing that the same process is applied to all the equations of the model. Consider b := P x ∈ Ω b and a test function where P p is the space of polynomials of degree at most p. Following a DG-FEM type recipe based on double integration by parts on each sub-domain (Dumbser and Facchini, 2016;Hesthaven and Warburton, 2007), we multiply the equation by ϕ and integrate over the domain to obtain the weak form: Double integrating by part the first derivative term permits to evaluate the coupling terms between the domains, following the discontinuous coupling method by (Eskilsson and Sherwin, 2006) where the boundary integral represents the flux through the boundary surface. The termP is a numerical representation of the boundary flux term and we chose to estimate it here with a centered flux. Thus, replacing the unknowns with a spectral/hp element approximation, the equation can be solved with just a projection matrix M , a derivative matrix operator D and a matrix to represent the coupling/border condition C. At discrete level the eq. (10) becomes For simplicity, we can define the matrix Q := D + C.
Combining the first derivative matrices, we discretize all the equations for the outer and the inner domain: and where D d is an extrapolation of the water or body elevation node by node. The acceleration is defined by Newton's second law in eq. (5). We denote with w the vector containing the Gauss-Lobatto-Legendre integration weights using the corresponding abscissas of the quadrature. The discrete formulation of eq. (5) thus read Provided that the matrix QD d M −1 Q is invertible, the acceleration can be pre-evaluated as: where we have defined the added mass M add

A Solitary Wave Passing a Fixed Pontoon
We here study the case of a weakly solitary wave propagating past a rectangular box (Lin, 2006). The case is as follows. Consider a pontoon of length L = 5m and draft T 0 = 0.4m in a flume of constant still water depth h 0 = 1.0m. The two wave gauges are located at G 1 = −31.5m and G 2 = 26.5m assuming the centre of the box located at x c = 0m. See figure 2. The incoming solitary wave (Bonneton et al., 2011) has a non-dimensional amplitude A/h 0 of 0.1. The simulation is done with a mesh of 25 elements in the free surface domain and 5 elements in the body domain. The total length of the numerical flume is 185m with the box in the centre of the flume. All computations are done with a polynomial order p = 3. Also visible in figure 2 (the dashed lines) are the extension of the transition layer. Figure 3 shows the solution at two different times -at the maximum and minimum wave elevation. We see that the flow is smooth over the domain boundaries and that there are virtually no oscillations present in the pressure acting on the body. The surface elevations registered by the wave gauges are presented in figure 4. In all, there is a good agreement between the RANS simulations and the Boussinesq model. The minor lag in time can partly be due to the fact that the initial solitary wave might differ between the models. The Boussinesq model slightly over-predicts the transmission and is not fully capable to capture the short reflected waves, but all the salient features of the solution are captured correctly. This model can be easily modified to solve the case of a fixed pontoon with different bottom shapes and see how this effect the transmitted and reflected wave. In particular we tested a triangular bottom, a round parabolic bottom and a flat bottom with a deeper submerged area. The non-squared pontoons are designed to have a submerged volume equal to the one of the first test in figure 2. The mesh is the same as for the previous pontoon test. We see that the shape greatly affects the profile of the inner pressure (figures 3f, 5c, 5f). However, the transmitted and reflected waves are very similar to the original test, as we can see from the gauges plot in figure 6a. For the fixed body setting the displaced volume (area in 2D) seems more important than the shape or the maximum depth. This is confirmed by the last test in figure 5g, where we have a flat bottom pontoon that reaches the same maximum depth of the triangular one thus having a larger displaced area. From the gauges in figure 6a we see a substantial difference as the reflected wave is higher and the transmitted smaller. The force applied vertically by the wave on the bottom of the body in figure  6b is very similar in all the tests with the biggest difference seen after the peak of the wave has passed the body.

Latching of a Heaving Box
As mentioned in the introduction: latching is a classical technique to increase the power production of a WEC by controlling the motion of the device (Budal et al., 1979). In short, the device is locked when the body velocity is zero and hold still for a period of time (latching period) and then released, resulting in a higher body velocity. Giorgi and Ringwood (2016) recently showed that the latching period obtained from linear models are different and sub-optimal compared to latching periods obtained from CFD simulations. The difference is mainly due to the fact that latching breaks the assumption of small amplitude motion. Hence, it is interesting to use the Boussinesq model for latching as there is no such assumption present.
Here we implement the latching algorithm as outlined in Giorgi and Ringwood (2016) and apply it to the case of a heaving box with a width of 6m and draft 5m in a water depth of 20m. The incoming waves can be regarded as linear (H = 2.23m) with a period of 8s, which is well outside the box's heave natural frequency of 5.5s. The results is seen in figure 7, where the latching have been activated at t = 40s. First we see that the latching, as expected, greatly increase the motion response of the box. The non-latched solution has a response amplitude operator (RAO) of roughly 1.2. The latched case, albeit the heave amplitude is not constant in time, clearly has a RAO above 2. We see some differences between the Boussinesq and the CFD results, mainly in phase, but again the Boussinesq model correctly captures the salient features of the solution. The latched response is not perfectly harmonic due to reflected waves and even some air entrainment (in the CFD case).

CONCLUSION
We have presented a nonlinear numerical model for wave-body interaction in shallow water. The model is based on depth-integrated Boussinesq-type equations, a computationally efficient method for wave propagation in near-shore regions. The model is based on the unified Boussinesq approach of Jiang (2001), in which the domain is decomposed into free surface and body domains. The equations are in all subdomains formulated in a depth-integrated manner and discretized using the spectral/hp element method. The result is a computationally efficient method. We have illustrated the use of the model for fixed and heaving truncated bodies. The pontoon test shows very good agreement between the Boussinesq model and the CFD solution by (Lin, 2006). We then used the Boussinesq model to investigate the influence of body shape on the transmission and reflection of a solitary wave, and found that the shape had small impact (for the same displaced volume). We also presented preliminary results regarding phase control (latching) in the Boussinesq framework. Compared to CFD computations there was a phase shift in the response. This issue is under investigation.
The continuation of this work aims to expand the framework to two spatial dimensions in the free surface plane and include more degree of freedom for the body. The final aspiration of this work is to bridge these fundamental developments to applications of engineering relevance.