Development and Validation of a Procedure for Numerical Vibration Analysis of an Oscillating Wave Surge Converter

During extreme sea states so called impact events can be observed on the wave energy converter Oyster. In small scale experimental tests these impact events cause high frequency signals in the measured load which decrease confidence in the data obtained. These loads depend on the structural dynamics of the model. Amplification of the loads can occur and is transferred through the structure from the point of impact to the load cell located in the foundation. Since the determination of design data and load cases for Wave Energy Converters originate from scale experiments, this lack of confidence has a direct effect on the development. Numerical vibration analysis is a valuable tool in the research of the structural load response of Oyster to impact events, but must take into account the effect of the surrounding water. This can be done efficiently by adding an added mass distribution, computed with a linearised potential boundary element method.Thispaperpresentsthedevelopmentandvalidationofanumericalprocedure,whichcouplesthe OpenSource boundary element code NEMOH with the Finite Element Analysis tool CodeAster. Numerical results of the natural frequencies and mode shapes of the structure under the influence of added mass due to specific structural modes are compared with experimental results.


Introduction
Wave Energy remains one of the largest unexploited renewable energy resources and has led to many conceptual developments of Wave Energy Converters (WECs) [1]. One such type and most advanced concept is the Oscillating Wave Energy Converter (OWSC) Oyster [2], developed by Aquamarine Power Limited (cf. Fig. 1). The Oyster OWSC consists of a buoyant hinged flap, which pitches back and forth in the near shore environment pumping water onshore to generate electricity in a conventional hydro power station [3]. It can be shown that diffraction governs the motion of Oyster, leading to a difference in water level on the seaward and landward faces of the flap. The resulting differential in hydrostatic pressure (cf. Fig. 2) creates an exciting torque around the hinge axis [4].
The operating location of the near shore offers several advantages; it is closer to land where the energy is generated and * Corresponding author. Tel.: +44 2890974012; fax: +44 2890974278.
E-mail address: p.schmitt@qub.ac.uk (P. Schmitt). needed; there is amplification of the horizontal particle motion and thus surge and pitch loads can be exploited; and there is a natural reduction in the extreme wave heights and thus extreme loads due to sea bed dissipation and depth limited wave breaking [5].

Loading
Survivability is a key driver in the development of WECs, they must be designed to be available, reliable and robust. Although the Oyster is able to duck underneath the incident waves and its position in the near shore environment lends itself to the natural filtering of extreme waves, the extreme significant wave height is still predicted to be 3 to 4 times the average [6]. The principal loads on an OWSC are described in six degrees of freedom (DOF), three forces Surge, Sway, Heave and three moments Roll, Pitch and Yaw as shown in Fig. 3. Pitch and Surge are the principal driving forces and the remaining seen as parasitic to the design. An accurate description of the loading is crucial. The loading can be split into two distinct categories:   loads, lower in magnitude but high in cycle count, approximately 100 million over a 25 year lifecycle of a WEC. These governing structural loads are due to the hydrodynamic pressure differential on the flap due to water surface elevation and can be described using experimental and numerical techniques with an understanding of the long term wave climate.
Extreme loads: These are the infrequent impulsive loads, nonlinear in nature and difficult to describe statistically and using experimental or numerical techniques. They are less understood in WECs compared to other marine structures which are typically designed to dissipate such loads.
An experimental programme by [7] to understand extreme loads observed the phenomenon of slamming as depicted in Fig. 4. Under certain conditions the flap can be seen to pitch landward under an approaching wave crest. Due to its buoyancy it then pitches back, the seaward face dries out in the wave trough and it accelerates towards the oncoming wave crest.
This fluid-body interaction then causes a local, short period high frequency impact. Structural dynamics cause the global foundation loads, which are measured by a six DOF strain gauged load cell, to be amplified and contaminated (cf. Fig. 5) [8].
Methodologies and techniques to describe and understand the load cases for design have been developed and many still rely on small scale experimental tank testing and numerical analysis. Since this experimental data is used for full scale design, high frequency loads and structural dynamics are of considerable importance to the structural integrity, durability and costs of a WEC. Design of scale models and the analysis methods have similitude to techniques of full scale design. The slamming phenomena has also been observed on the full scale prototype device of Oyster installed in EMEC Orkney. Thus understanding and describing the local flap loads, global foundation loads and the structural response of the system to loads is vital.
Vibrational analysis at both prototype and experimental scales can be used to increase the understanding of the structural dynamics of Oyster and thus the definition of extreme loads. Once natural frequencies are understood it enables the scaled data to be either filtered, to remove the contamination, or, if deemed a full scale phenomenon, considered in the design of the structure or eliminated. From the images in Fig. 4 it can be seen that the main challenge is the fluid-structure interaction and the dynamic influence of the surrounding water. It is possible   to couple structural finite element simulations with fully nonlinear viscous computational fluid dynamics solvers [9], but the computational burden is still extremely high for practical engineering applications, like the design process of a WEC. [10] presented a simpler method for lattice structures like off-shore platforms, where the effect of the water was included in a modal analysis by using Morrison's equation to compute an added mass for each beam element. Unfortunately, the OWC cannot be described well using Morrison's Equation since its mode of operation is dominated by diffraction, as described in detail in [4]. Many finite element codes allow to solve the Euler equation for closed domains, like pipes filled with a fluid or beams in a filled tank [11]. These methods usually require the discretisation of the entire fluid domain [12] and sometimes cannot take into account free surface conditions. Boundary element codes that efficiently solve a linearised, potential flow problem have been used in marine engineering for decades, but are still a useful tool [13] and actively developed [14]. These tools allow to derive added mass and damping with very little computational burden. Especially for vibration related problems, with typically relatively small motion amplitudes, the linearisations employed are mostly valid and this class of tools is thus well suited for coupling with structural finite element simulations.
At the time of writing no method is known to the authors which couples the added mass due to the fluid and wave interaction of a OWSC with standard finite element analysis (FEA) to perform a modal analysis of a OWSC structure.

Objectives
This paper describes the development and validation of a numerical tool that couples two OpenSource software tools, Code_Aster, an industry developed FEA tool and NEMOH, a Boundary Element Method (BEM) code. As a first test we will calculate the (water level dependent) added mass for the device for different rates of immersion and secondly carry out a numerical vibration analysis to understand natural frequencies and mode shapes. A simplified version of an Oyster scale model (cf. Fig. 6) is used to obtain experimental data for comparison, as described in Section 6. This model is unable to pitch about the pivot axis in order to keep the experimental and numerical implementation simple and limit sources of errors.

Added mass
Any moving or deforming body in contact with fluids experiences loads caused by these [15]. Intensive research to understand and determine these loads is found in marine engineering under the term of Added Mass. Added mass occurs as soon as an acceleration a is imposed to a fluid and the kinetic energy E kin related to the fluid motion changes [16]. Physically the added mass can either be visualised as the necessary work which has to be done to achieve this change of E kin or as an additional mass attached to the body in contact with the fluid [16]. Assuming a body with a fluid-structure interface (wetted surface) Ω moves through the surrounding fluid with a constant velocity v Eq. (1) describes the kinetic energy related to the fluid motion with ρ being the constant density of an incompressible fluid [15].
Due to the constance of the fluid-structure interface the integral  Ω dΩ can be approximated with an invariant value N [16] so that Eq. (1) can be simplified to Assuming an accelerated body with a time variant velocity described as d v /d t, the kinetic energy also becomes time variant ( d E kin/d t) which leads to an additional drag force F drag experienced by the moving body specified by As ρ · N describes a mass, the characterisation of added mass as an additional mass accelerated with the body is straightforward [15]. It should be noted that the term added mass does not necessarily include only mass (m added ) values but also embraces moments of inertia (I added ) and masses multiplied by lengths [17]. Therefore the term Added Mass Coefficients 11 · · · · · · · · · · · · · · · · · · m 22 · · · · · · · · · · · · · · · · · · m 33 · · · · · · · · · · · · · · · · · · I 11 · · · · · · · · · · · · · · · · · · I 22 · · · · · · · · · · · · · · · · · · I 33 is used in the following to describe a tensor, made up of the added mass, inertia and cross coupling terms (omitted here) acting in the three translational and rotational degrees of freedom.
For some three dimensional geometries such as ship hulls analytical solutions for the added mass coefficients are available [16]. With a stepwise calculation considering linear potential flow firstly CM is determined for the two dimensional cross section of the investigated geometry and secondly the third dimension is considered by the so called J factor (three dimensional correction factor) which is assumed to be constant along the hull [16].
The calculation of added mass for complex three dimensional geometries requires numerical methods. The so called surface panel method determines added mass from fluid forces on the structure [16]. Again considering potential flow theory, the velocity potential ϕ can be defined for each panel of the fluid-structure interface Ω so that with the help of the Bernoulli equation pressure p is defined as Integration over all surface panels then leads to the fluid related force F stated in Eq. (6) [12].
With a known panel acceleration the derivation of CM is straightforward. If the free surface of the investigated system includes the appearance of fluid waves those have to be considered in the free surface boundary condition [12]. Defining the wave amplitude A, the fluid depth h and the surface elevation η = h + A the free surface equation of motion dη /dt = − dA /dt leads to the surface boundary condition for the potential flow dA dt = ∂ϕ ∂η (7) which directly affects the derivation of CM [12].

Natural frequency
The natural frequency describes the free vibration frequency of a mechanical system imposed by an initial excitation impulse [18]. Applying Newton's second law of motion one obtains the governing equation of motion for the derivation of the natural frequency of a system This formulation can be substituted into Eq. (8) and, by neglecting damping, leads to the Eigenvalue problem with the characteristic equation

Eigenmodes of a body in a dense fluid
In Eq. (10) the square of the natural frequency ω represents the Eigenvalue with the corresponding Eigenvector {Ψ } which represents the mode shape of the system [19]. Neglecting additional hydrostatic stiffness the structural properties of a mechanical system have to be modified when considering the surrounding fluid of a structure by adding the added mass coefficient matrix to Eq. (10). Hence a direct influence on the natural frequency by CM can be observed (cf. Eq. (12)).

Modal analysis
The term Modal Analysis embraces the derivation of the three characteristics: natural frequency, mode shapes and structural damping factors [19]. The modal analysis can either be performed theoretically or experimentally.
The aim of the experimental modal analysis is the determination of the frequency response function (FRF) describing the relationship between the vibration response of the investigated system (X (ω)) and the excitation (P(ω)) by defining a function dependent on the exciting frequency (H(ω)) [19].
The FRF measures the response rate in terms of displacement, velocity, or acceleration of a structure to an excitation force [20].
To determine frequency response functions, data for the excitation and the vibration response are required. The two most common experimental setups are either shaker or impact hammer tests [19,20]. During shaker tests a fixed excitation location has to be chosen and vibration responses are measured at varying positions. From the response function the natural frequencies of the structure can be derived. For impact hammer tests fixed response measuring locations (e.g. fixed accelerometer positions) are used to capture the vibration response to spatially distributed excitations [19]. To perform impact hammer tests, the structure being investigated is hit with a hammer and allowed to vibrate freely in its natural frequency. Fig. 7 shows an acceleration time trace obtained from an impact hammer test, as recorded in the experiments described later in Section 6.
By choosing specific impact locations one can excite the structure in a desired degree of freedom and mode shape. From the measured response the damped natural frequency is obtained. The damped natural frequency is a function of undamped natural frequency and damping, though often assumed to be identical.  Fig. 6. It should be noted, that the grey PVC plates on the front were not considered in the simulations, since their effect is negligible. Different parts, represented by element colours, were meshed separately as second order volume elements. Almost all connection sets are based on forcing surfaces on a slave face to deform identical to the touching nodes of the coupled 3D elements. Earlier numerical investigations showed, that the dynamic properties are dominated by the weight of the flap structure and the flexibility of the load-cell and baseplate. The   Table 1. The accuracy of the structural model is later assessed in Section 7.1 by comparing results for the numerical and experimental modal analysis without the effect of surrounding water.

Nemoh BEM
The OpenSource BEM solver NEMOH performs calculations on a two dimensional mesh describing the shell of the structure being investigated. An example of a boundary mesh describing a simple OWSC model is shown in Fig. 9. With user defined boundary conditions such as water depth, motion frequency and position of rotation axis hydrodynamic coefficients are calculated. Among others (e.g. wave induced excitation forces) the added mass coefficients are calculated and provided in matrix form (cf. (14)). This calculation of the coefficients is based on the pressure distribution of the investigated structure (cf. Section 4.1). Assuming a rigid body, one gets added mass coefficients matrices for the six different degrees of freedom surge, sway, heave, roll, pitch and yaw. The rigid body assumption enables one to calculate relevant added mass coefficients for the vibration analysis considering the first modes for the rotational DOF's. However, to get reasonable results for modes of higher order, complex mode shapes need to be considered for the added mass calculation with the BEM code. The current work focuses on these  basic modes, but the method will be extended to arbitrary mode shapes in the future. When calculating added mass coefficients for the 40th scale static OWSC FEA model it is assumed that the flap has the biggest effect on the added mass. Hence added mass coefficients are only calculated for the flap shell to keep computational effort low.
Added mass coefficients plotted over a frequency bandwidth in-between 0 and 16 Hz are shown in Figs. 10 and 11 for a half and a fully submerged flap, respectively. The typical form of the added mass over frequency can be observed in all figures. For very low frequencies, the added mass is extremely high, decreases sharply to a minimum for slightly increased frequencies before slowly increasing towards a constant value at high frequencies. Some simulations show physically unreasonable results at low frequencies, but for the current application the region below 8 Hz is not of interest.
For the application of an added mass matrix to the structural simulation it has to be verified that the natural frequencies are high enough to be almost constant for the frequency range of interest. The added mass matrices calculated by the BEM with the boundary conditions listed in Table 2 for a half respectively fully submerged flap shell are shown in matrix (14) and (15).

Coupling methodology Derivation of added mass from radiation pressures
The radiation pressure distribution given by NEMOH for each DOF is used to derive a discrete mass element distribution.
NEMOH internally defines surface elevation η as η = −A · sin(ωt) (16) where A represents the unit wave amplitude and ω the wave frequency [14]. The coupling between surface elevation and pressure on the body surface p leads to Eq. (17).
WithP as absolute pressure (per unit velocity and wave frequency) and phase angle Θ. Due to the fact that the unit velocity differs for translational and rotational motions added mass has to be calculated in two different ways for translational and rotational modes. For translational motion the use of the NEMOH output pressures is straightforward (cf. Eq. (18)) since a unit velocity of 1 m /s is considered to normalise pressure. For rotational motion the unit velocity of 1 rad /s is used for the normalisation, hence the given pressureP has to be divided by the distance to the rotation axis r (cf. Eq. (19)).
These two pressures can then be used to calculate a pressure related force vector ⃗ F (cf. Eq. (20)) and the local added mass vector ⃗ m added (cf. Eq. (21)) by multiplying with the surface normal vector ⃗ n (∥⃗ n∥ = Φ with Φ as surface area) and dividing by the incident wave frequency ω, respectively.
From the local added mass vector the local added inertia ⃗ I added can then be determined by multiplying with the square of the distance to the rotation axis (cf. Eq. (22)).
To do so a convention for the distances to the rotation axis r has to be declared. Values for r are defined as With ⃗ c as vector of the centre node coordinate with its components in surge (⃗ c surge ), sway (⃗ c sway ) and heave (⃗ c heave ) as shown in Fig. 12.

Automated mass derivation
In order to implement a reasonably high resolved added mass distribution to the FEA model a methodology is needed which uses the NEMOH output data as an input and gives out added mass vectors which can be defined as discrete (mass) elements by Code_Aster. Such a methodology has been developed by using several python functions coupled with the Code_Aster environment. In the first place the NEMOH output file containing the coordinates of each pressure point namely the centre of each panel ⃗ c surge , ⃗ c sway , ⃗ c heave (cf. Fig. 12), the absolute pressure per unit velocityP and the phase angle Θ has to be converted into a python compatible array. To derive ⃗

Code_Aster interface
Code_Aster provides an added mass element type. This element type implements directionally dependent mass properties and is assigned to mesh nodes. A distribution of these discrete mass elements, are passed to Code_Aster as python lists. The projection of this data onto the 3D ''modal mesh'' is performed in python. The ''modal mesh'' here means the three dimensional mesh of the structure on which the modal analysis (or any other Finite Element Analysis) is performed.

Experimental tank tests
The natural frequencies for the three DOF's roll, pitch and yaw of a static OWSC are obtained experimentally by choosing three different impact locations. Acceleration and strain gauge data are recorded with two accelerometers and a six degree of freedom load cell, respectively (cf. Fig. 13). As a test case the flap has been tested at different water levels. Water level W1 in Fig. 13 represents a submerged load cell, W2 in Fig. 13 represents a half submerged upright flap and W3 a fully submerged upright flap model. Table 3 shows the natural frequencies for vibrations in roll, pitch and yaw for the different water levels. The considerable effect added mass has on the structural dynamics can be observed by

Table 3
Resulting natural frequencies normalised with frequencies in a dry environment.

Roll
Pitch Yaw comparing the obtained Eigenfrequencies at varying water levels. The fully submerged flap has an Eigenfrequency of only 33% of the dry structure. It should also be noted, that the decrease of natural frequency of the system is not simply proportional to the increase in water level. Raising the water level from half to fully submerged does not lead to a change in frequency of 100% but varies for each mode.

Initial structural model
To test the numerical model on physical reliability it is common practice to compare numerical and physical natural frequencies [19,20]. Since the natural frequency is affected by the structural dynamic properties this measurement can show discrepancies between the two (cf. Section 4.1). The structural behaviour of a physical as well as numerical model is mainly affected by the material properties and connections sets between different parts. In the numerical model the results are also affected by the discretisation of the structure thus the mesh properties (number of mesh nodes, number of mesh elements, etc.). To obtain the natural frequency of the numerical model a FEA modal analysis has been carried out. With the help of experimental modal analysis of the dry model the FEA model of Oyster was validated. Results for the natural frequency in a dry environment are presented in Table 4. The maximum error is 5% for the yaw mode, which is much less than the effect of the added mass. The model is thus deemed

Preliminary investigations
The mass of the surrounding water which has to be accelerated when the structure moves can simply be described as a body attached to the actual structure having certain mass and inertia properties for the six different modes surge, sway, heave, roll, pitch and yaw. Since the first Eigenmodes of this specific model consist of rotation of the quasi undeformed flap structure around the loadcell, it might be believed, that the implementation of a single lumped mass with the mass and inertia properties of the surrounding water into the FEA structure would yield good results. Tests were performed using the added mass coefficients from matrix (14) and (15). The resulting frequencies for a fully and half submerged flap are shown in Table 5. Results show that the simple assumption of single lumped mass does not capture the effect of the surrounding fluid. Deviation between experimental and numerical results is very high. Possibly the coupling of a significant mass and inertia in a single node of the flap structure is not suitable. Since the model should be used Table 5 Relative deviation between experimental and numerical results for the eigenfrequency obtained from simulations with a single lumped mass and varying water levels.

Roll
Pitch Yaw in the future for arbitrary structures and arbitrary modeshapes, where the entire distribution is required, this simplification was not investigated further.

Verification of the added mass distribution
Due to the failure of a single lumped mass implementation to mimic added mass effects of an OWSC the methodology presented in Section 4.1 is used. Simulations for a half and a fully submerged static 40th scale flap model are performed. The pressure distribution over the flap shell for the three DOF's roll, pitch and yaw are calculated with NEMOH. Although not directly used in the FE Analysis, added mass is also calculated for the three translational motions surge, sway, heave as well as added inertia for the three translational motions roll, pitch and yaw. To gain confidence in the mapping algorithm and implementation in the software, these values are compared to the NEMOH outputs (cf .  Tables 6 and 7   and coloured according to computed pitch and roll modes. The normalised frequencies can be seen in Fig. 16.

Discussion
The mass and inertia values for the DOF's surge, sway, heave, pitch and yaw show good agreement between the original NEMOH data and the transformed added mass/inertia applied to the structure. The maximum deviation is only 6% (cf. Tables 6 and 7). Natural frequencies obtained for these modes with the numerical modal analysis agree very well to experimental data, with maximum deviations of about 7%.
Discrepancies of 15% and 50% for the added inertia in NEMOH and on the structural model are found for the half and fully submerged roll mode, respectively. The resulting natural frequencies for a rolling flap show differences of 8% and 22% for the half and the fully submerged flap (cf. Fig. 16). Since the considerable error in the application of added mass does not lead to an error of the same magnitude in the prediction of the Eigenfrequency, it could be argued that in this particular mode structural stiffness and mass dominate the vibration. Hydrostatic stiffness, which is zero for the yaw mode and also small in pitch, could be considerable for roll motion. At the same time the added mass is also much larger for pitch (1.9 kg m 2 ) and yaw (0.92 kg m 2 ) modes than for roll (0.24 kg m 2 ), when the flap only accelerates water through the displacement of the relatively small bottom face of the flap. The natural frequencies from the numerical simulation is also expected to be slightly higher than the damped natural frequency obtained from the experiments. Table 6 Results for added mass and added inertia for a half submerged flap from NEMOH and Python/Code_Aster: Dominating Added Mass Vector Components highlighted Bold.

Conclusion
This paper presented an experimental and numerical investigation into the effect a surrounding fluid has on a vibrating OWSC structure. A method to couple an efficient BEM potential flow solver with a structural finite element solver was presented and validated with experimental data. The performed verification study shows good results and proves the functionality of the developed methodology.
The main conclusions are: • The surrounding fluid must be considered in any dynamic structural analysis of an OWSC.
• A single lumped added mass is not able to mimic the physical reality.
• Coupling linearised, potential BEM solvers with FEA tools is an efficient way to accurately model the effect of added mass.
For the presented verification study only data for an upright static flap model at two different water levels has been used. The method should be extended in the future to take into account hydrostatic stiffness, which would also enable the investigation of moving flaps and probably improve accuracy for the roll mode.