Introducing non-rigid body structural dynamics to WEC-Sim

This paper describes the development of a structural dynamics add-on to WEC-Sim, an open-source code dedicated to the dynamic analysis of Wave Energy Converters (WECs). When calculating the dynamic response of a body, WEC-Sim by default uses a rigid body dynamics approach. Such an approach ignores the potential effects of structural deformation on the WEC, which can in turn affect e.g. the distributed loads across the WEC and / or the individual (point) load sources that depend on the dynamic response of the WEC. Following a similar approach to tools used in the offshore wind industry, a structural dynamic add-on was developed using Code_Aster as the Finite Element (FE) solver to enable coupled hydro-elastic, timedomain analysis. The add-on was developed and tested using an example Oscillating Wave Surge Converter (OWSC) WEC model, currently being developed as part of the H2020 MegaRoller project. In the examples studied, the inclusion of structural dynamics is shown to affect the estimated peak Power Take-Off (PTO) loads, with variations in PTO force of over 10% being observed when structural dynamics are considered in the analysis. Keywords— Wave Energy Converter (WEC), load analysis, structural analysis, coupled analysis, WEC-Sim.


A. Background
The Horizon 2020 MegaRoller project aims to design, build and validate a high performance, cost-efficient and reliable 1MW-scale Power Take-Off (PTO) that can be integrated into Oscillating Wave Surge Converter (OWSC) designs. In this context, the knowledge of localised effects related to both the environmental conditions and the WEC response is of crucial importance to the design, implementation and integration activities to be completed. One of the consortium's key focuses is therefore to increase the understanding of the wave-structure interaction problem via the development of a fully-coupled, nonlinear loads model, suitable for the assessment of distributed loads and for performance and structural analysis from an early design stage up to the transition to detail design. In this context, this paper follows several previous studies which have investigated the wave-structure interaction of OWSC's. For example [1] reviews some of the work completed in the period 2012 -2016.
WEC-Sim [2] is a time-domain tool dedicated to the prediction of loads for Wave Energy Converters (WEC's). In its original formulation, and when calculating the dynamic response of a body, WEC-Sim takes a rigid body dynamics approach. In such an approach, each hydrodynamic body is treated as a single point with mass and inertia properties, and the loading outputs from a WEC-Sim assessment are restricted to point loads acting on the entire device and a limited set of distributed pressures (namely the hydrostatic and Froude-Krylov components). Without a more detailed representation of the WEC structure it is not possible to calculate loads at specific locations and / or at key structural features (e.g. stiffeners, spars).
In addition, the current approach neglects the potential effects of structural deformations on the WEC dynamics, which can potentially affect e.g. the distributed loads across the WEC and / or the individual (point) load sources that depend on the dynamic response of the WEC. As an extreme example, hydrodynamic loads may excite fundamental frequencies of the structure, leading to resonance and potentially the self-destruction of the WEC.
In order to accurately enable the consideration of structural deformations of the prime mover of a WEC, such as the main body (flap) of an OWSC, a structural dynamics add-on is required.

B. Overview of the proposed approach
Previous work by Y. Guo et al. [3] has introduced a form of structural dynamics to WEC-Sim. A generalised modes approach was used, where the system's degrees-offreedom (DOFs) are reduced based on the structure's modal response. However, this requires prior knowledge or calculation of the WECs' mode shapes and limits the structural flexibility to a limited number of DoFs.
The offshore wind industry provides an example of where loads prediction tools that include structural dynamics have been developed. Typically, a coupled, hydro-elastic, time-domain approach is followed; see e.g. HydroDyn [4]. Such approach often uses a Finite-Element (FE) method to solve the structural dynamics.
The influence of the coupling on the system's response was evidenced in previous studies. Coupling should in this context be understood as the consideration of structural dynamics in loads calculations. For example, in [5] the loads on offshore jacket foundations were found to vary by up to 25% when coupled analysis is considered, compared to those derived using an uncoupled approach.
Following a similar methodology, this paper describes the development of a structural dynamic add-on using a FE solver in conjunction with WEC-Sim, enabling hydroelastic, time-domain analysis of WECs. Compared with generalised modes approach, such methodology allows the modelling of more generic structures where the mode shapes may not be known in advance.
The development of the add-on first considered the use of the FE model as a post-processing tool (for both static or dynamic analyses). Ultimately, the dynamic FE model was coupled to the WEC-Sim model in a similar fashion to other tools (e.g. HydroDyn).
The paper is organised in seven main sections. Following this introduction, Sections II and III describe the overall WEC-Sim and FE models developed, respectively. Sections IV, V and VI then detail three case studies of increasing complexity: Case studies 1 and 2 use an uncoupled, post-processing approach to the test the FE model setup, firstly using a stationary WEC with a static FE solver (Case study 1, Section IV), an then a dynamic model capable of predicting the time-dependent response of the WEC (Case study 2, Section V). Finally, the dynamic model is coupled to the WEC-Sim loads calculations, enabling coupled hydro-elastic analysis (Case study 3, Section VI). The paper is concluded in Section VII with a shortlist of recommended next steps.

A. Overall description of the WEC-Sim model
The MegaRoller WEC (see Fig. 1) can be described as an OWSC with a modular PTO solution, where hydraulic piston pumps have an interface with the prime mover (flap) via a twin drivetrain located at each end of the flap. The device operates in nearshore regions at depths of approximately 8 to 20 metres. It is anchored to the seabed and, depending on tidal range, it is mostly or fully submerged during operation. A series of devices can be deployed in an array to create a wave farm. Since the WEC is constructed as a modular individual unit, there is no technical upper limit to the number of devices that can be used in an array.
A WEC-Sim model of the OWSC was developed as part of the MegaRoller project to enable loads calculations, which will in turn inform the design of the PTO units. At a high level, the numerical model consists of a hydrodynamic body representing the gravity base, connected to the seabed via a fixed constraint, and to the flap via a pitching constraint (mimicking the bearings), as illustrated in Fig. 2. Each PTO subsystem includes a translational PTO block. The model is described in detail in [6].
In this paper two design situations were considered. In a first case study (Case study 1, see Section IV), the bearing constraint (see Fig. 2) is locked, restraining the flap in all degrees of freedom (DoF), simulating a fault scenario leading to the locking of the PTO units and / or bearings. The model is run for an example irregular sea state, defined as a JONSWAP spectrum with significant wave height of 2.75m and a peak period of 11.5s. In a second step (Case studies 2 and 3, see Sections V and VI), the bearing constraint is unlocked and an alternative fault condition is considered whereby the PTO unit on the left side of the flap (PTO1) features a damping coefficient ( 1 ) five times higher than the PTO unit on the right side (PTO2). The model is run in the same sea state as in case study 1 but with the waves approaching at a 20-degree angle to the normal of the flap. The Case studies differ in the approach taken, with Case study 3 using a coupled approach whereas Case study 2 uses an uncoupled, postprocessing approach.
The key properties of both the models and the input environmental conditions are summarised in Table I, for each of the case studies presented in Sections IV to VI.

B. Post-processing of WEC-Sim outputs
A structural dynamics tool requires distributed pressures (and/or loads) as an input to capture any difference in loading across the flap. Previous work by Michelen et al. [7] has demonstrated a method for extending the default WEC-Sim outputs to include distributed diffraction and added-mass pressures, using the output from the NEMOH BEM code [8]. Extending the regular wave case presented in [7] to generic, irregular sea state applications, the diffraction pressure (1) and addedmass pressure (2) can be calculated as following: where is the diffraction pressure on cell , , and are the frequency, direction and phase of the ℎ wave component, . is the time and is the diffraction pressure coefficient for cell , output from NEMOH. And, where is the pressure associated with the addedmass on cell , ̈ is the acceleration for DoF number and k,i is the radiation pressure coefficient for cell and DoF number , output from NEMOH.
When using the convolution integral formulation, WEC-Sim uses the added-mass coefficient at infinite frequency ( ∞ ). As this coefficient is calculated by NEMOH only for the whole body (and not for each cell), the added-mass pressure was approximated using the distributed added-mass coefficient at the highest frequency considered (i.e., = 4 rad/s in (2)). Similarly, the method for calculating radiation damping pressures presented in [7] is not applicable when using the convolution integral formulation. Distributed viscous drag forces were also unavailable. As an initial approximation, the radiation damping and viscous drag distributions were therefore applied as uniform pressure fields , with the magnitude calculated using the WEC-Sim point load for these force components divided by the surface area of the flap. Future development work could focus on developing a distributed formulation for these two force contributions -see also Section VII. Fig. 3 illustrates the diffraction and added-mass pressure distribution for Case study 2 (see Table I). The 20deg offset of the incoming wave direction causes an asymmetric distribution of the diffraction pressure along the width of the flap. Overall, diffraction and added-mass pressures vary in magnitude across the surface of the flap. This may lead to non-uniform deformation and unequal loading, which may lead to effects that can only be captured when using a non-rigid body structural dynamics approach.

A. Structural model
A structural dynamics tool requires the development of a separate FE model of the MegaRoller device, which can use the distributed pressure loads (see also Section II.B.) as an input and calculate the distributed displacements and nodal forces. Code_Aster [9], a widely used open-source code capable of a range of different types of analysis (e.g., static, dynamic, modal, nonlinear), was used as the FE tool for this assessment.  As a first approximation, the MegaRoller flap was approximated to a structure discretised into a number of steel panels and internal stiffeners. All panels were assumed to feature a thickness of 25mm. The density of individual sections was then adjusted to match the overall mass, centre of gravity and inertia as listed in Table I. The panels are thin compared to their length and width, so Code_Aster's 'DKT' shell elements were used for the structural mesh. A relatively coarse mesh was applied to enable solution times comparable with those of WEC-Sim. It should be noted that the FE model was conceived with the goal of predicting displacements and nodal forces, not detailed stresses, and so a coarse mesh was considered appropriate for this stage. The mesh used, defined with an average element size of 3m 2 , is illustrated in Fig. 4. The baseline FE model was assigned linear-elastic steel material properties (Youngs Modulus E = 210GPa, Poisson Ratio v = 0.3). The overall effect of different materials in the structural response was then investigated by varying the Youngs Modulus during the case studies presented in Section IV to VI.
The following subsections provide a high-level description of the static and dynamic approaches that were considered in the FE model. The level of coupling is also introduced, where the FE model can be used following uncoupled (post-processing only) or coupled approaches.

B. Static approach
In a static analysis the inertial and structural damping effects are neglected and Hooke's law (4) applies: where is the position, the stiffness matrix of the structure, and is the resultant force.
The static approach can either be used to analyse single instant in time, or to approximate a time-series using a series of static 'snapshots'. By removing the timedependent terms from the finite-element equations, this approach allows large transients to be reduced to a much smaller number of analysis steps, reducing computational effort and/or allowing for more detailed structural analysis. However, it should be noted that the approach assumes that the excitation frequency is much lower than the natural frequencies of the structure (typically a maximum ratio of 1/3). For WECs, which are often designed to operate near-resonance, such assumption may be a limiting factor.
As velocity and acceleration are not estimated, the static method can only be applied in a post-processing approach, where the structural model is solved separately to WEC-Sim, as illustrated in Fig. 5. This enables the calculation of distributed loads and deformations but does not allow the deformations to influence the loads calculations. Such approach is similar to the method demonstrated in [10] but using WEC-Sim as the hydrodynamic solver instead of WaveDyn.

C. Dynamic approach
In a dynamic FE analysis, the equation of motion (5) is solved to estimate the full transient response, including both inertial, structural damping and structural stiffness effects. This is equivalent to the approach used by a rigidbody solver such as WEC-Sim, except that for rigid-bodies the structural stiffness and damping matrices are zero.
where ,,̈ are the position, velocity and acceleration of the body, is the mass matrix, the structural damping matrix, the structural stiffness matrix, and ( ) is the (dynamic) resultant force. Using the position, velocity and acceleration estimates, the dynamic method can be coupled to WEC-Sim. In a coupled approach between the hydrodynamics and the structural model, distributed pressure loads are derived from WEC-Sim and applied to the FE model at each time step. The FE model is then used to estimate the position, velocity and acceleration at the flap centre of gravity and PTO interfaces. These are then passed back to the WEC-Sim model to calculate the loads at the next time step. Such a sequential method can be termed "loose coupling", following the same approach used to couple MoorDyn to WEC-Sim [11] and as described in more detail by Ahamad et al. [12]. The procedure is illustrated in Fig. 6. An alternative approach, the sensitivity to which could be  This coupling allows structural deformation of the flap to influence the loads calculation, in particular uneven deformation across the flap which can result in unequal PTO loads. However, the method still relies on rigid-body hydrodynamic coefficients and cannot, for example, update the diffraction force coefficients based on the deformed shape of the flap. Future work could implement a more complex methodology using e.g. a database of hydrodynamic coefficients (see also Section VII).
The main disadvantage of the dynamic method is the increased solution time due to the need for small timesteps, both to ensure stability and to adequately capture the transient response.

IV.
CASE STUDY 1 -STATIC POST-PROCESSING ANALYSIS As a first step towards implementing a structural dynamics module in WEC-Sim, the first case study focused on an uncoupled, static analysis of the MegaRoller OWSC. All DoFs were restrained at the bearings and PTO interfaces (see Fig. 4) to simulate the flap being locked in position.
The 'STAT_NON_LINE' solver in Code_Aster was used to implement a static analysis. The model was run using both the baseline steel material properties and a 'flexible' case, where the material Youngs Modulus was reduced to 5GPa to assess and highlight any resulting structural deformation.
As discussed in Section III.B, a static FE solution is not time-dependant. As such, it does not require the analysis of every time instant from the WEC-Sim, so long as a sufficient number of events are selected to adequately capture the structural response. In this case study, an iterative method was applied to downsample a set of structural time instants based on the applied surge force, limiting the maximum difference between the full WEC-Sim solution and the linearly interpolated structural timepoints to 250kN (~2.5% of the peak force). Fig. 7 shows the final selection of 458 samples, compared to a total of 50,000 steps in the WEC-Sim time signal. Fig. 8 compares the linearly interpolated reaction force in surge between the structural model (in blue) and the WEC-Sim (in orange) output. The two signals are generally well matched, verifying the accuracy of the distributed diffraction, Froude-Krylov and hydrostatic pressure loads. The distributed loads in the structural model tend to slightly over-estimate the WEC-Sim surge load at the peak events. This is likely to be due to the differences in the mesh used between the two models. Future studies could aim to investigate the sensitivity of the model to the mesh density.
The structural results provide fully-distributed nodal forces and displacements, enabling more detailed loads assessments at specific points of the structure (when compared to the rigid-body solver). Fig. 9 illustrates the displacement in surge at the top of the flap, estimated for a steel structure (in red) and for a more flexible (generic) material ( = 5 , in blue). Using steel properties, the flap structure is relatively stiff, with a maximum surge displacement at the top of the flap of 0.06m (see scale on the right-hand-side of the figure). However, and as expected. using a more flexible material increases the magnitude of this displacement, with a maximum surge displacement at the top of the flap of 1.18m (see scale on the left-hand-side of the figure). Fig. 10 illustrates the deformed shape for the flap constructed using the flexible material. In addition to the global bending of the flap, deformation of each panel between stiffeners is also visible. This local deformation is  The second case study adapted the structural model to include a dynamic solver, enabling the prediction of the flap's angular motion, velocity and acceleration. The FE model boundary conditions were modified to allow pitch rotation at the bearings and apply point loads at the PTO interfaces.
The Code_Aster dynamics solver 'DYNA_NON_LINE' was used to model the transient response -see (5 In this study an iterative approach was used to determine damping coefficients such that high frequency oscillations were damped, and the model was stable, resulting in coefficients = 0.01 and = 0.0025. Future work could seek to refine these estimates based on a more detailed assessment of geometry and material propertiessee also Section VII. The same post-processing approach as used in Case study 1 was taken (see Fig. 5), where distributed pressure loads and point loads at the PTO interfaces were calculated after completion of the WEC-Sim simulation and then applied to a separate structural model. The case study was run for the 'steel' and 'flexible' ( = 5 ) materials. As an initial verification exercise of the model setup and implementation of the dynamic solver, Fig. 11 compares the flap rotation estimates between the structural model (in black) and the WEC-Sim model (in blue), for the baseline steel material. The result is generally a close match to the WEC-Sim rigid-body solution, with small differences mostly associated with the model setup (e.g., the relatively coarse structural mesh, small differences in mass and inertia properties).
The potential benefits of considering structural dynamics is highlighted when analysing the velocity at the PTO units. Fig. 12 compares the velocity estimates from the left PTO (i.e. the one with the highest damping coefficient 1 ) between the structural model using steel (in blue), the structural model using a more flexible material (in red) and the rigid WEC-Sim model (in black). For the stiff, steel material, the predicted velocity matches closely the WEC-Sim rigid-body solution. When a flexible material is used, the asymmetric loading due to the 20-degree offset in the wave direction and the unequal PTO damping (see also  Table I) leads to a significant bending of the flap, as illustrated in Fig. 13, and results in noticeable differences in the velocity at the PTO interface. Although the material stiffness used in this example is artificially low and unlikely to represent an actual WEC, such results demonstrate the potential impact of distributed loads on the structural response. Case study 2 assessed the capability of the FE model to estimate distributed position, velocity and acceleration results across the flap. However, the differences in motion at key interfaces due to the structural deformations were not accounted for in the load's calculations in WEC-Sim. In the final case study, the 'Flap' hydrodynamic body in WEC-Sim (see Fig. 2) was customised to enable the coupling with the FE structural model.
To assess the influence of the coupling on the MegaRoller WEC's response, the coupled WEC-Sim model was first run using steel material properties. Fig. 14 compares the rotation of the flap and the PTO response between the coupled approach (Case study 3, in blue), the post-processing (Case study 2, in dashed blue) and the WEC-Sim rigid-body solution (in black) outputs.
In this example, the coupling between the hydrodynamic and structural models leads to a larger range of PTO forces and flap rotation. This is likely due to the increased local velocity caused by structural deformation at the PTO interface. The same coupled model was used to conduct an initial investigation into the sensitivity of the PTO to the flap's material stiffness. A less stiff material with Young's Modulus E = 100GPa was used to represent a hypothetical composite material (e.g. GFRP). Fig. 15 compares the PTO force and velocity for the PTO unit with the highest damping coefficient ( 1 ) between the coupled approach for the stiff (steel, in blue) and flexible ( = 100 , in ref) materials and the WEC-Sim (in black) outputs. It should be noted that, for the original WEC-Sim results, the PTO velocity output is the same for both PTO units as the flap remains rigid throughout the simulation.
The inclusion of structural dynamics has a noticeable effect on the PTO response, as shown in Table II in terms of peak and RMS PTO loads: the absolute peak PTO load over the 100s simulation increased from 713kN and 792kN (approximatively 11%) when comparing the coupled steel and uncoupled cases. When comparing the RMS of the PTO force signal, an increase from 232kN to 258kN, again approximately 11%, is observed.
In addition, the use of a more flexible material lead to an overall decrease of the PTO loads, with the peak PTO load reducing by approximately 12% (from 792kN to 700kN) between a stiff material a more flexible material. Further investigations during the MegaRoller project will aim to test this sensitivity for a wider range of load cases (see also Section VII).  Finally, the coupling of WEC-Sim to a structural model, solved at every time step of the simulation, leads to a significant increase of the computational time (c.12,310s in this example), approximately 13 times longer when compared to a nonlinear WEC-Sim simulation (c. 920s), and more than 300 times longer when compared to a linear WEC-Sim simulation (40s). Future work will seek to optimise computational efficiency in order to increase the suitability of the tool for use in analysing multiple load cases (see also Section VII). However, even in its current form the analysis of selected design load cases is feasible and solution time is still significantly faster than, for example, typical CFD simulations.

VII. CONCLUSIONS AND NEXT STEPS
This paper documents the development and testing of a novel structural dynamics add-on for WEC-Sim, which enables the consideration of structural deformations. Three steps of increasing complexity were considered and illustrated using the MegaRoller OWSC as a case study, demonstrating the potential impact of a coupled hydroelastic model on the WEC's response.
Firstly, Code_Aster was used to develop a static FE model of the MegaRoller WEC. The model was used in a post-processing approach, where distributed pressures were calculated in WEC-Sim and applied separately to the structural model. Reaction forces were shown to closely match the WEC-Sim solution, verifying the distributed pressure estimates. The preliminary model enabled the estimation of distributed WEC loads and displacements (whereas the rigid-body WEC-Sim solution does not predict any local displacements).
Secondly the structural post-processing module was updated to consider a dynamic analysis. Comparisons of the flap rotation time-series output showed an overall good agreement between the rigid-body WEC-Sim solution and the structural model outputs. The dynamic analysis enabled the estimation of the local deformation and dynamics, specifically at the PTO interfaces.
In the final case study, a coupled hydro-elastic model was developed, where the results from the dynamic FE model are fed back in WEC-Sim's loads calculations at each time step. For the example load case studied, the inclusion of coupled structural dynamics was shown to change the predicted maximum peak PTO loads by up to 11% compared to the WEC-Sim output.
It should also be noted that the case studies presented in this paper used a passive control strategy with a constant damping coefficient for the PTO units. More complex cases, for example including both PTO stiffness and damping terms or step changes in the response, could further exacerbate the structural response. In addition, the radiation and viscous drag force components were modelled as uniform pressures as a first approximation. In reality, the distribution of these forces across the flap could further affect the structural deformation. As an immediate next step, it is envisaged a formulation distributing the radiation and viscous drag pressures will be developed, replacing the uniform loads applied in this study. Further developments could aim at the following: • The level of hydro-elastic coupling could be increased to extend the suitability of the tool for other WEC devices, for example by updating WEC-Sim's nonlinear hydrostatic calculations at each time step or by implementing a database of hydrodynamic coefficients which are updated based on the deformed shape at each time step. • The damping coefficients considered in the Rayleigh approximation could be refined based on a more detailed geometry and material properties. • The sensitivity of the results of this study to the mesh density and 'loose' coupling algorithm used could be assessed. • Finally, the computational efficiency of the tool could be increased by better integrating Code_Aster with WEC-Sim and providing a more automated workflow.
In the context of the MegaRoller project, the new tool will be used in DLC assessments to predict loads on the twin drive-train PTO design, assess different PTO fault scenarios and further investigate the sensitivity of the WEC to material stiffness.