Two-Dimensional Modelling of a Quayside Floating System

This paper studies the behaviour of a quayside floater oscillating in front of a vertical dike. In order to study the floater motion and the impact of the dike on the floater, a linear analytical model based on 2D potential flow theory in intermediate water depth conditions and a numerical model resolving 2D Navier–Stokes equations are developed. Physical tests performed for different floater dimensions in a wave tank are used as references for the analytical and numerical models. The comparison of the results obtained analytically, numerically and experimentally leads to the validity domain of the potential model. A correction of this model is proposed, based on the optimization of the radiated coefficients, and a quadratic drag term is added according to Morison equation. The impact of the different parameters of the system on floater behaviour is considered. Results show that the draft has the most important impact on floater motion.


Introduction
Production of energy is a major environmental challenge. The development of renewable marine energy systems has been increasing during the last years. Many technologies of Wave Energy Converters (WEC) have been developed and may be firstly classified in offshore, near-shore and on-shore (quayside) systems [1,2]. Although wave energy potential is lower at coast than offshore, quayside systems offer several advantages. Indeed, although near-shore energy is about 10% lower than offshore energy, the quayside systems also have the advantage of protecting the dikes against oceanic storms by limiting the waves overtopping [3,4]. In addition to this advantage, costs in terms of maintenance and grid connections are largely reduced compared to offshore systems and waves become more regular when reaching the shore. Recent studies show a trend to use near-shore or quayside systems in order to reduce maintenance and construction costs [5,6]. Moreover, the presence of the dike can be used to improve the energy produced thanks to the wave's reflection. The present study is situated in this context.
If motion of floating buoys and WEC efficiency are often obtained from experimental and complete numerical approaches, preliminary design is generally done by simplified numerical or analytical models [7]. Indeed calculation time of numerical modelling can be long for preliminary design [8,9]. Linear models can then offer an alternative approach.
The model described by Andersen and Wuzhou [10] is based on Green's functions and eigenfunction in order to reduce the calculation time of the previous models. Some years later, simplified formulas for added mass and damping coefficients were proposed by Drimer et al. [11].
Results were confirmed by comparing an analytical model to a linear numerical one. Berggren and Johansson [12] also proposed calculation of the radiated coefficients of a two-body system (composed of floating cylinder connected to a submerged plate) by using the matched eigenfunction expansions method. They obtained a reasonable agreement between the radiated coefficients of their system and those of a single floating cylinder. This study has been used by Drobyshevski [13] for the calculation of zero and infinite added mass. Results were compared to the linear numerical model results, and it was concluded that this model can be used to find an approximation of damping and added mass. Based on these previous results, Zheng et al. [14,15] and Zhang and Zhou [16] proposed to calculate hydrodynamic coefficients for offshore and quayside point absorber based on linear wave theory. These studies are essentially two-dimensional. The development of codes based on Boundary Elements Methods (BEM) enables the computation of hydrodynamic coefficients in 3D configurations [17][18][19].
Linearised wave theory has been used to study the interaction of floating structures with a vertical wall in various configurations. Elchahal et al. [5] studied the effects of reflection coefficient on a floating breakwater performances using a linear model. They highlighted the important role of the toe clearance, draft and floater's beam length in the transmission coefficient for large toe clearance. Bhattacharjee and Soares [20,21] studied the wave amplitude and the reflection coefficient between a vertical wall and a rectangular floater for a step type bottom using linear wave theory. Forces on the floating structure and the wall are studied as a function of the different parameters for toe clearance at least three times the draft. In particular, toe clearance and floater beam length can be modified to optimize free surface elevation between the dike and the floater, and wave loads on the wall.
In order to improve the validity of analytical models, non-linear effects can be added. Following this idea, Bonfiglio et al. [22] and Henning [23] evaluated numerically added mass and damping, including viscous effects. This method gives an improvement of results based on traditional evaluation of radiated coefficients. Confrontation of numerical and analytical models has been carried out by Bhinder et al. [24], who proposed a comparison of a potential time domain viscous model with different numerical models. The hydrodynamic coefficients are calculated by AQUAPLUS code [25], and compared with a viscous analytical model. They have proven that adding a Morison like quadratic drag coefficient (calibrated from the floater motion) is an improvement of the model. Corrections of added mass depending on the draft of the floater have been proposed by Koo and Kim [26] for an offshore point absorber. It appears from this model that corrected added mass in low frequencies differs from the numerical values due to strong bottom effect.
However for quayside systems, adding a corrected term may be insufficient. The aim of the present work was to study the behaviour of a two-dimensional quayside floating buoy in intermediate water depth with three different approaches: analytical, numerical and physical modelling. Comparison of the three models is used to determine the validity domain of the analytical model. The problem is simplified by considering a rectangular buoy, moving only along a vertical wall breakwater. The impact of the dike on the floater behaviour is also studied. Innovating points of this study are the comparison of the three models, the proximity of the system with the dike and the intermediate water depth condition.
The results of present work can be used for quayside floating systems to be established on intermediate water sites with important exposition to waves, such as the site of Esquibien in Brittany (France) where a 300 m long and 10 m high dike is present [27].
The model based on linear wave theory is described in Section 1. Numerical model (which solves the 2D Navier-Stokes equations) and experimental tests are considered in Sections 2 and 3, respectively. The results and the three models' comparison are given in Section 4. The last section is devoted to the conclusion.

Floater Motion
The system is composed of a bidimensional floating buoy oscillating in heave motion under monochromatic waves action, along a vertical wall. The floater is a homogeneous rectangle of mass m, width 2l, height H f and draft Γ, installed at a toe clearance D from the wall. The water depth h is uniform along the flume. The origin of the coordinate system is at the free surface level at rest, half-width of the floater, z coordinate is measured positive in the upward direction and x coordinate positive in the incident waves direction. The domain is virtually separated in three zones: zone (1) is situated between the dike and the floater, zone (2) is under the floater and zone (3) is the wave inlet zone (see Figure 1). The floater motion along the z-axis is expressed as follows: with Z(t) the vertical position of the floater's centre of mass, g is the gravity acceleration and S w the wetted surface of the floater [28]. Using the linearised potential flow theory at first order, pressure is obtained with Bernoulli's equation: where Φ is is the potential from which the velocity can be driven. The equation of motion thus becomes: where Π(t) represents the hydrostatic, i.e., buoyancy, forces. In the present configuration, gravity and hydrostatic forces are expressed as a spring force with a stiffness K Ar To express the hydrodynamic forces, the Haskind decomposition [29,30] of the potential flow is used. The potential Φ(x, z, t) is thus expressed as a sum of diffracted and radiated components (Equation (5)): where Φ I represents the incident potential, Φ D the diffracted potential (corresponding to the waves reflected by a fixed floater) and Φ R represents the radiated potential (which corresponds to the waves generated by a moving floater in still water). The hydrodynamic forces are then expressed as a sum of the excitation force F EX , computed from Φ D and Φ I (Equation (6)) and the radiation force F R , computed from Φ R (Equation (7)). Waves are assumed regular and monochromatic Considering monochromatic waves forcing, the variables are expressed in the frequency domain where the radiation force amplitude is expressed as the sum of an inertial and damping term: with c m the added mass, c a the damping coefficient andZ the complex floater motion amplitude. According to Equation (3), the floater heave motion can finally be written as follows: withF EX the complex amplitude of the excitation force. The real floater motion amplitude is then defined as A f = |Z|, (|Z| is the modulus of the complex numberZ).
To solve Equation (9), the hydrodynamic coefficients c m , c a and the excitation force amplitude |F EX | must be determined. To identify these coefficients, the values of the incident, diffracted and radiated potentials have to be assessed, according to Equations (6) and (7). Different identification methods of the hydrodynamic coefficients have been tested by Yang et al. [31] and Andersen and Wuzhou [10]. In the present study, the method of Zheng et al. [15], with eigenfunction matching, is used and described in the following subsection.

Eigenfunction Matching Method
The diffracted and radiated potentials are written as a sum of N terms composed of unknown coefficients (sum of fundamental and evanescent modes) to be determined [14,15]. N is the number of terms necessary for the desired precision. These coefficients are calculated so as to verify the boundary conditions and the continuity equations between the different zones. The diffracted and radiated depend on the considered zone (1, 2 or 3) and on the floater motion mode (heave, roll or switch). Here, only the heave motion is considered. Thus, the coefficients A 1n , A 2n , A 3n and B 2n represent the coefficients to be determined to obtain the diffracted wave potential, while the coefficients A 1n , A 2n , A 3n and B 2n represent the coefficients to be determined to obtain the radiated potential. finally, The expressions of different potentials in the different zones of the flow will also be expressed according to the eigenvalues λ n and β n given by: λ n tan(λ n h) = −ω 2 /g for n = 2, 3, ...., N

Expressions of the Radiated Wave Potential
-Zone 1: -Zone 2: -Zone 3: The boundaries and continuity conditions are used to solve the unknown coefficients. Those conditions take into account the dike, i.e., a total reflection of the wave. Continuity conditions of pressure and velocity are described between zones 1-2 and 2-3. A detailed description of this model can be found in Zheng et al. [15].
To find the unknown coefficients of the diffracted and radiated potential, each part of the equation obtained from the boundary and continuity condition is multiplied by its eigenfunction and integrated over its respective depth interval, at the float boundaries x = ±l. Only the N first terms of the infinite series are used, to deduce a linear system of 4N equations (Equations (19) and (20)), where X D and X R are the solution vectors for diffracted and radiated coefficients, respectively, S, G D and G R are the matrices obtained with boundaries and continuity conditions. The solution of the linear system given by Equations (19) and (20) is coded under PYTHON using the Generalized Minimal RESidual (GMRES) algorithm with the Incomplete Lower Upper factorisation with 0 fill-in (ILU(0)) preconditioner [32]. The nonlinear dispersion relations given by Equations (10) and (11) are solved by Newton-Raphson algorithm.
where the components of X D and X R are The code is validated by calculating the hydrodynamic coefficients for a case showed by Zheng et al. [15], i.e., h/l = 6, Γ/h = 0.4, D/h = 0.2 and N = 30 ( Figure 2). The convergence has been quantified by varying the number of terms N in the series. The relative difference between 10 and 20 terms is 1.8% and between 20 and 30 terms is 0.81%. For all calculations, a minimum of 30 terms is used.

Numerical Model
The problem is solved numerically using ANSYS Fluent software. The complete 2D Navier-Stokes equations are solved for a two phases flow (air + water) using the volume of fluid (VOF) method. The turbulence is solved using the standard k − model. The solution is calculated with the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithm due to Patankar [33]. The floater motion, wave amplitude and hydrodynamics coefficients are measured for each configurations, for 15 to 20 waves cycles depending on the wave period. The calculation of the incident and reflected waves are made using Goda and Suzuki [34] method with two numerical probes.

Convergence Study
Geometry and mesh are obtained with the ANSYS Workbench software. The mesh is chosen rectangular for all the tests, and is defined for the three zones defined in Figure 1. A mesh convergence has been studied following the convergence procedure of Celik et al. [35]. Three meshes types have been tested (Table 1). For each mesh, the size of the mesh is half the size of the previous tested mesh. The Courant-Friedrichs-Lewy condition (CFL) is defined as a characteristic velocity multiplied by the time step and divided by the cell dimension along the x axis. Values of CFL are calculated for zone 2 with the floater velocity amplitude and for zone 3 with the wave phase velocity ω/k, where ω = 3.74 rad.s −1 and k = 1.79 m −1 .
For each mesh, the floater motion and the wave amplitude are measured. Using the convergence procedure of Celik et al. [35], the relative error is about 0.91% for floater motion and 1.1% for wave amplitude between Mesh 1 and Mesh 2. The convergence index for Mesh 1 is about 1.14% for floater motion and 2.18% for wave amplitude. Finally, the mesh convergence study allowed us to observe that it is the Mesh 2 which offers the good compromise between calculation time and accuracy. Consequently, the size of this mesh will be adopted for all the simulations which follow. The time step is fixed at 0.0025 s.

Hydrodynamic Coefficients Calculation
The calculation of the numerical radiated coefficients has been carried out according to Henning [23] and Cozijn et al. [36]. The floater velocity is imposed with a User Defined Function (UDF) from ANSYS Fluent, following Equation (21): where A f f is the forced floater's amplitude. Calculated radiation force F z (t) = |F z | sin(ωt + φ) acting on the floater are decomposed in Fourier series. Hydrodynamic coefficients are then calculated by identifying Fourier coefficients a 1 and b 1 , to the motion equation terms Equation (9) in the absence of body mass. Therefore, the added mass (c m,num ) and damping (c a,num ) coefficients can be obtained numerically by: Figure 3 presents the comparison between the radiation force directly extracted from ANSYS Fluent software and the one deduced from the numerical coefficients c m,num and c a,num identified above by the expressions (22) and (23). It should be noted that these expressions were obtained for the case of floater 5 (Table 2) and for a period of wave T = 1.29 s.
A good agreement between the numerical and the fitted radiation force with a relative error less than 2% for all simulations is found.

Experimental Conditions
The The incident and reflected wave amplitudes are calculated using the method of Mansard and Funke [37] with four resistive probes placed in the flume. The maximum uncertainty of wave amplitude measurements is 1 mm, corresponding to a maximum uncertainty of 2% for all tests. A 2000 × 2000 pixels CCD camera Basler aca2040 is placed in front of the floater to record its motion using the software HIRIS R&D Vision at a frequency of 20 Hz. The floater is painted in black and a white target is added to follow its trajectory. The maximum error made on the floater motion is smaller than 1 mm including MATLAB post-processing, corresponding to a maximum uncertainty of 3% for all tests. The erro bar (barely visible) is plotted in Figure 8c for only one test (the last one) since the error is the same for all tests.
A parametric study is performed to highlit the influence of the floater's dimensions and hydrodynamics forcing on the floater behavior. Thus, the floater's draft, height, width and its distance to the dike have been varied, according to Table 2 Table 2). A total of 61 tests has been carried out. Five independent dimensionless numbers are defined: Γ/l represents the floater aspect ratio, Γ/h the water depth effect, kD the effect of the toe clearance on the floater behaviour, kh the wave dispersion and A i /h the wave non linearity, with A i the incident wave amplitude and k = 2π/λ where λ is the wavelength. The dimensionless numbers Γ/l, Γ/h , kD, kh and A i /h are summed-up in Table 2.

Impact of the Dike
In order to study the influence of the dike on the floater dynamics, the motion amplitude of floaters 1 and 5 is calculated analytically with and without the dike (as performed by [14,15]) for the same wave conditions ( Figure 5). These results are compared to numerical data with the dike. Firstly, the proximity of the dike modifies the hydrodynamic coefficients leading to a decrease of the natural frequency. In presence of the dike, analytical values of the hydrodynamic coefficients are about one and a half times larger than the ones without the dike at small frequencies as shown in Figure 6. For higher frequencies, the effects of the dike on the floater disappear, except for the added mass coefficient, the damping and diffraction coefficients have close values with or without the dike when ω Γ/g ≥ 1. However, radiation coefficients calculated with the numerical model are larger than analytical ones, leading to a strong decrease of the natural frequency for both floaters and a decrease of the floater motion amplitude for floater 5. Thus, the non-linear phenomena, modelled in numerical simulations (but not in analytical simulations), strongly increase hydrodynamic coefficients.
Furthermore, analytical results with the dike present a second resonance at a dimensionless frequency ω Γ/g close to one. This resonance corresponds to a piston mode of the fluid in the gap between the floater and the dike. The value of this peak frequency is in good agreement with the estimation of Molin [38] for the piston mode in moonpool. Figure 7 represents the motion amplitude of floater 5 for a toe clearance varying from D = 0.01 m to 0.10 m. Analytical results are compared to experimental measurements. When the toe clearance is increased, the natural frequency of the floater motion is slightly decreased analytically although no clear difference is observed on experimental results. In the same time, the floater motion amplitude is decreased.  Moreover, for large values of the frequency, the amplitude of the floater motion, measured experimentally, strongly decreases. This decrease is enhanced as the toe clearance is increased. Thus, the floater's operating range of frequencies is reduced for large toe clearance. Therefore, the natural frequency for a quayside system is decreased compared to an open sea system. This decrease is enhanced when taking into account non-linear phenomena. Moreover the motion amplitude is decreased when increasing the toe clearance.

Impact of the Water Depth and Draft
At first, the amplitude of motion of the floater 1 is calculated for different water depths with a constant aspect ratio D/Γ = 0.2. Figure 8 presents the comparison between the results of the numerical, analytical and physical models of the dimensionless floater amplitude A f /A i . Figure 8 shows a satisfactory agreement between the experimental and numerical results. Except for the largest water depth, the natural frequency is similar for both models. However, the amplitude of motion is slightly underestimated in the numerical model. The analytical model fails at predicting the natural frequency of the floater with an error of about 17%. On the other hand, the water depth has very few impacts on the floater behaviour and no impact on the accuracy of the analytical model. Dimensionless floaters motion with different drafts is calculated with the different models ( Figure 9). As for floater 1, a good agreement is found for floater 5 between physical and numerical models ( Figure 9d); thus, only numerical simulations have been performed from floaters 2 to 4. When the draft is increased from Γ = 0.05 m to Γ = 0.25 m, the natural frequency is decreased from 6.97 rad.s −1 to 3.74 rad.s −1 . Moreover, the resonance amplitude increases for increasing values of the draft. For drafts varying from Γ = 0.05 m to Γ = 0.20 m, i.e, D/Γ ∈ [0.05, 0.2], numerical dimensionless amplitudes present a secondary peak at a frequency lower than the natural frequency. For Γ = 0.25 m, i.e., D/Γ = 0.04, those two peaks seem to merge. At last, for increasing values of the draft, the difference between the numerical and analytical models increases. Both motion amplitude at the resonance and the natural frequency of the floater in the linear model, strongly increase when increasing the draft. Thus, as the draft is increased, dissipation due to non-linear phenomena increase. To understand the difference between linear analytical and numerical models, the spanwise component of vorticity is calculated for floaters 1 and 5 ( Figure 10) at the same frequency and D = 0.01 m. Vorticity field shows strong vortices generated at the edges of the floater dues to the non-linear effects in the floater movement. As the draft increases, the intensity of the vortices increases from Ω max = 5.2 s −1 for Γ = 0.05 m (Figure 10, left) to Ω max = 16.5 s −1 for Γ = 0.25 m (Figure 10, right). The vortices shed at the edges of the floater are stronger for larger drafts leading to stronger non-linear effects and energy dissipation. Moreover, as the draft is increased, a combination of a Couette and Poiseuille flow within the gap is set-up leading to strong dissipation. The viscous force thus increases with the aspect ratio Γ/D. This may explain the differences between the inviscid linear analytical model and numerical simulations.

Impact of the Incident WAVE amplitude
The floater motion has been measured for different wave amplitudes. Only experimental values are considered in this part. The response of the floater A f /A i with ω/ω 0 is shown in Figure 11 Present results show that the water depth h has no significant impact on the floater motion amplitude whereas the draft Γ and the toe clearance D have important impacts. Depending on the system dimensions, it has been observed that the analytical model does not represent properly the floater behaviour for all tested cases. Indeed, the resonance frequency and amplitude of the floater motion may be overestimated by the analytical model, compared to experimental and numerical results.

Conclusions
The behaviour of a floater in heave motion in front of a vertical dike has been studied with three approaches: analytical, numerical and experimental. The analytical model is based on the solution of the heave potentials following Zheng et al. [15] and using the Haskind decomposition. Amplitude of the wave potential modes are solved using the eigenfunction matching method. Hydrodynamic coefficients are then calculated from heave potentials. The results of the analytical model have been compared with numerical simulations and experimental tests. The 2D Navier-Stokes equations are solved using the ANSYS Fluent software. Numerical results are in good agreement with experimental tests performed in a wave flume. The influence of the water depth, the draft, the toe clearance and the wave amplitude on the floater motion is studied.
The proximity of a dike modifies the natural frequency of the floater. Moreover, resonances of fluid in the gap between the floater and the dike are generated at dimensionless frequencies close one. As the draft is increased the natural frequency of the floater get closer to the piston mode frequency. However, the linear model fails at predicting both the natural frequency and resonance amplitude of the floater for small values of D/Γ. Indeed, non-linear dissipation due to vortex shedding and viscous terms increase as D/Γ is decreased. Vortices generated at the edges of the floater are then stronger and viscous dissipation in the gap increases.

Funding:
The project is co-financed by CEREMA, and the Normandy council and European Union through ERDF funds (NEPTUNE project).