Flexible multibody dynamics modelling of point-absorber wave energy converters

wave loadings


Introduction
With the exhaustion of fossil fuel resources and growing demand for energy consumption, many countries (e.g. the UK and China) are making considerable efforts to exploit renewable energy resources [1e3].As an inexhaustible and environmentally-friendly energy resource, ocean wave power, which is extracted from ocean waves through WECs (wave energy converters), is highly valued by coastal countries [4,5].The worldwide resource of wave energy is estimated to be in the order of 1 TW [6].For European waters, the wave resource is estimated to be capable of covering more than 50% of the total power consumption [7].
According to the ways used to capture the energy of the waves, WECs can be roughly categorised into eight groups [8], i.e. attenuator, point absorber, oscillating wave surge converter, oscillating water column, overtopping device, submerged pressure differential, bulge wave and rotating mass.Fig. 1 presents the schematics of the eight types of WECs.Attenuator (see Fig. 1a), which operates parallel to the wave direction and rides the waves, captures energy from the relative motion of the two arms as the wave passes them.Point absorber (see Fig. 1b), which absorbs energy from all directions through its movements near the water surface, converts the heave motion of the floater into electrical power.Oscillating wave surge converter (see Fig. 1c), the arm of which oscillates as a pendulum mounted on a pivoted joint in response to the movement of water in the waves, extracts energy from wave surges and the movement of water particles within them.Oscillating water column (see Fig. 1d), which uses a large volume of moving water as a piston in a cylinder to compress and decompress the air column, turns a turbine at the top of the column to generate electricity by the movement of the air.Overtopping device (see Fig. 1e) captures the sea water of incident waves in a reservoir above the sea level, and then releases the water back to the sea through a conventional low-head turbine which generates power.Submerged pressure differential (see Fig. 1f), which is typically located near shore and attached to the seabed, uses the alternating pressure in the device caused by the motion of the wave to pump fluid through a system to generate electricity.Bulge wave (see Fig. 1g) device consists of rubber tubes filled with water and moored to the seabed.The water enters through the stern and the passing wave causes pressure variations along the tube, creating a 'bulge'.The bulge grows as it travels through the tube, gathering energy which can be used to drive a turbine at the end of the tube.Rotating mass (see Fig. 1h), which uses two forms of rotation to capture energy by the movement of the device heaving and swaying in the waves, drives an eccentric weight to rotate to generate electricity.Among these types of WECs, the point-absorber WEC has the advantage of ease of both fabrication and installation, and therefore has been widely studied [9,10].
Due to their relatively small size, the amount of energy that the conventional single point-absorber WECs can capture is relatively small when compared to other types of WECs.This limitation can be overcome by using multiple point absorbers, which consists of several floaters.Using several floaters also contributes to increase system reliability and stabilise the platform.
The key component of a WEC is its PTO (power take-off) system, which converts the wave induced motion of the floaters into an electrical power output to the grid.A reliable structural model, which is capable of accurately determining structural dynamic responses under different load cases and limit states, is crucial in the structural design of the PTO system.It can also provide valuable information regarding the locations of hot spots which should be reinforced and potentially monitored through an integrated structural health monitoring system [11].
For the structural model, the FEA (finite element analysis) model has been widely used for solving complex static and dynamic problems encountered in engineering, due to its high fidelity and flexibility.It has also been extensively applied to the structural modelling of renewable energy devices.Authors have applied FEA to the fluid structure interaction of both HAWT (horizontal axis wind turbine) blades [12] and VAWT (vertical axis wind turbine) rotors [13] as well as structural optimisation of VAWT composite blades [14].The application of FEA to marine structures can be found in Refs.[15e17].
As the PTO of a WEC is a complex system comprising multiple bodies, accurate modelling of PTO also requires MBD (multibody dynamics) analysis.The combination of FEA and MBD is the socalled FMBD (flexibility multibody dynamics), in which the interconnected bodies of the multibody system are considered flexible (i.e. each body is not rigid and allowed to be deformed).
This paper develops a new analytical model particularly relevant to WECs and presents in detail the development of a high fidelity MBD-FEA that can be applied not only to WECs but marine energy systems of a similar class.The floater arm tip displacement and velocity obtained from the FMBD model are validated against the results obtained from the analytical model of the structural modelling of the NOTC 10 kW multi-point-absorber WEC.After the validation, the FMBD model is further used to calculate the stress distributions, fatigue life, deformations, modal frequencies and modal shapes of the structure.The two approaches presented in this paper can be of value to the further development of marine energy technologies.
This paper is structured as follows.Section 2 presents the NOTC 10 kW multiple-point-absorber WEC.Section 3 presents the analytical model of the point-absorber WEC.Section 4 presents the FMBD model.Results and discussion are provided in Section 5, followed by conclusions in Section 6.

NTOC 10 kW multi-point-absorber WEC
The WEC used in this study is the NOTC 10 kW multi-pointabsorber WEC (see Fig. 2), which is designed by NOTC (National Ocean Technology Centre) and manufactured by THOECL (Tianjing Haijin Ocean Engineering Corporation Limited).It follows the principle of multiple point absorbers, utilising a ship-type platform for the prototype.It comprises six cone-shape floaters, which can heave in the wave and extract the power at a total capacity of 10 kW.With the help of floater arms, the motion of each floater can be transferred to the piston, pumping a hydraulic cylinder to output high pressure hydraulic oil.The hydraulic oil is transported to the hydraulic motor through pipelines, and then the hydraulic motor outputs energy to the generator to produce electricity.The diameter of the floater and the length of the floater arm are 1.52 m and 2.40 m, respectively.Fig. 3 depicts the 3D (three-dimensional) geometry model.Fig. 3a presents the 3D geometry model of the NTOC 10 kW multipoint-absorber WEC, and Fig. 3b shows a close view of the floater.As can be seen from Fig. 3a, the WEC comprises six identical floaters.Therefore, only a single floater is needed in the modelling.As can be seen from Fig. 3b, with the help of the floater arm, the motion of the floater is transferred to the piston, pumping the oil in the hydraulic cylinder.The hydraulic cylinder is connected to the bracket, the bottom surface of which is fixed to the platform.The floater arm, piston, hydraulic cylinder and bracket constitute an MBD (multibody dynamics) system, the assembly geometry model of which is presented in Fig. 4.

Analytical model of point-absorber WECs
An analytical model of point-absorber WECs is developed in this work, and the schematic of the analytical model is presented in Fig. 5.As can be seen from Fig. 5, the floater arm is represented by Bar AD having an angle of q with respect to the Z axis.Point A is a pivot point offsetting from the water surface with a height of h.Point D offsets from the water surface with a height of z and is hinged to the floater.The PTO, which consists of piston, hydraulic cylinder and hydraulic oil, is represented by Bar BC having an angle of 4 with respect to the Z axis.Point B is a pivot point offsetting from Point A with a distance of a, and Point C is hinged to the floater arm.The b 1 , b and c denote the lengths of AC, AD and BC, respectively.
According to the geometrical relations given in Fig. 5, the length c and floater top location zcan be respectively expressed as: The wave load calculation and kinetic model of point-absorber WECs are presented below.

Wave load calculations
Considering the heave motion of the floater, the perturbation velocity potential field of the floater in the wave satisfies the following conditions: where f is the disturbed potential of the water; subscript in denotes the potential of the incident wave; n is the unit normal direction of the floater surface; u 3 is the moving speed of the floater in a vertical direction (along Z-axis in Fig. 5); n 3 is the component of n in a vertical direction; t is time; g is the gravity acceleration; z is the vertical position of any water mass point.
For monochromatic incident waves, the velocity potential of the wave is given by: where a is the wave amplitude of the waves; Refgdenotes the real part of a complex value; F 0 in is the frequency-domain wave potential of monochromatic incident waves with unit amplitude; u is the circular frequency of monochromatic incident waves.
In case of infinite water depth, F 0 in in Eq. ( 6) can be expressed as: where k is the wave number.
r in Eq. ( 7) is given by: where c is the direction angle of incident waves; x andy are coordinates.
The wave surface equation is given by: where x is the lift of free surface.
When the floater is in the equilibrium state, the frequency response of the floater is identical to the wave frequency.Thus, u 3 in Eq. ( 4) can be expressed as: where U 3 is the amplitude of velocity u 3 .Due to its linearility, the perturbation velocity potential can be decomposed to: where F 0 is the potential of the diffracted wave in the frequency domain corresponding to the incident waves with unit amplitude; F 3 is the potential of the radiation wave in the frequency domain in a vertical direction when U 3 ¼ 1.
F 0 and F 3 in Eq. ( 11) satisfy the following conditions: v vn The linear hydrodynamic pressure at an arbitrary point on the floater surface can be expressed as: Substituting Eqs. ( 6) and ( 11) into Eq.( 18) gives: The linear wave loads in the direction of floater heave, F 3 , can be expressed as: where SW denotes the wet surface of the floater.Substituting Eq. ( 19) into Eq.( 20) yields: Define the force coefficient F 0 3 , added mass l 33 and damping coefficient m 33 as: where Imfgdenotes the imaginary part of a complex value.With the help of Eqs. ( 22), ( 23) and ( 24), the heaving wave loads F 3 in Eq. ( 21) can be rewritten as:

Kinetic model of point-absorber WECs
The total kinetic energy of the wave device can be expressed as: where E 1 is the kinetic energy due to the rotation of the floater arm (Bar AD) around Point A, E 2 is the kinetic energy due to the rotation of the hydraulic cylinder and piston (Bar BC) around Point B, E 3 is the kinetic energy of the line motion of the hydraulic cylinder and piston, E 4 is the kinetic energy due to the heave motion of the floater.E 1 , E 2 , E 3 and E 4 in Eq. ( 26) are respectively given by: where J is the moment of inertia of Bar AD around point A;q is the angle between Bar AD and the horizontal direction;J P is the moment of inertia of PTO (Bar BC) around point B;4 is the angle between PTO BC and horizontal direction;m P is the mass of PTO;c is the length of PTO between points B and C in Fig. 5;m is the mass of the floater; z is the vertical position of the floater.
According to the law of conservation of energy, the following equation can be obtained: where r w is the density of water; A W is the water-plane area of the floater;F 0 2 is the frequency-domain wave exciting force action on the floater;l 22 is the added mass of the floater;m 22 is the wavemaking damp of the floater;K P is the spring stiffness of the PTO;m P is the damp of the PTO.
ÀrA W ðz À z 0 Þ in Eq. ( 31) is the static restoring force caused by the change of draught of the floater, and its value is derived from the difference between the gravity and buoyancy of the floater.K P ðc À c 0 Þ is the linear restoring force during the movement of the PTO.m P _ c is the linear damping of the PTO, and the instantaneous output power of the wave energy converter is given by: With the help of Eqs. ( 26) to (30), Eq. ( 31) can be rewritten as: Recalling Eqs. ( 1) and ( 2), the following equations can be derived through differentiation: Introducing: Eqs. ( 34) to (37) can be then rewritten as: Considering the mass of the PTO is far less than the mass of moving parts gives: Substituting Eqs.(42) to (47) into Eq.(33) yields:

FMBD (flexible multibody dynamics) model of pointabsorber WECs
A FMBD model of point-absorber WECs is developed using the transient structural module of ANSYS software package, a widely used commercial finite element software package.The adjacent bodies of the multibody system are connected through joints, such as revolute and translational joints.All the bodies are considered flexible and the stresses within each body are evaluated using FEA.The geometry, material, mesh, joints and load boundary conditions used in the FMBD model are presented below.

Geometry
The floater arm, piston, hydraulic cylinder and bracket constitute a multibody dynamics system, the geometry model of which is presented in Fig. 6.

Material
The floater arm, piston, hydraulic cylinder and bracket of the NOTC 10 kW prototype are fabricated with different types of steel; however, details of the materials used in the fabrication are not allowed to be disclosed in this paper for confidentiality reasons.Therefore, in this study, the floater arm, piston, hydraulic cylinder and bracket are assumed to be made of S355, which is a widely used steel for marine structures.The mechanical properties of S355 steel are listed in Table 1.

Mesh
The types of mesh used for FEA can be roughly categorised into two groups, i.e. structured and unstructured mesh.Compared to unstructured mesh, structured mesh has the advantage of providing more reliable results.Additionally, filling the same volume of a structure, the number of elements required by structured mesh is much fewer than those required by unstructured mesh, and therefore structured mesh can save much computational time when compared to the unstructured mesh.For this reason, structured mesh is used in this study.In order to determine the proper mesh size, mesh independence/convergence study needs to be carried out [19,20].In this case, a point load of 1 kN is applied to the floater arm tip.Four mesh sizes are investigated, i.e. 0.040 m, 0.020 m and 0.010 m, and 0.005 m.The associated total number of elements and the calculated maximum von Mises stresses are presented in Fig. 7 and Table 2.As can be seen from Fig. 7 and Table 2, the maximum von Mises stress converges at a mesh size of 0.010 m, with a relative difference (2%) when compared to further mesh refinement with a mesh size of 0.005 m.Therefore, 0.010 m is deemed as the appropriate mesh size.The created mesh is depicted in Fig. 8.

Joints
In the multibody dynamic system, the relative motion between two adjacent bodies is constrained by joints.In this study, three types of joint are used, i.e. fixed, revolute and translational.The fixed joint constrains all DOFs (degrees of freedom).The revolute joint constrains three relative displacement DOFs and two relative rotational DOFs, leaving only one rotational DOF available.The translational joint constrains the two relative displacement DOFs and three rotational DOFs, leaving only one displacement DOF available.

Table 1
Material properties of S355 [18].The joints used in this study are presented in Fig. 9 and Table 3.
As can be seen from Fig. 9 and Table 3, the bottom surface of the bracket is fixed by using a fixed joint.A revolute joint is used to connect the bracket and the hydraulic cylinder, allowing relative revolution between them.A translation joint is defined between the hydraulic cylinder and the piston, allowing the piston to move along the hydraulic cylinder.A revolute joint is used to connect the piston and floater arm, allowing relative revolution between them.A revolute joint is defined at the end of the floater arm, allowing the floater arm to revolve around the centre of its end cylinders.
Additionally, a spring joint, which is characterised by PTO stiffness and damping, is defined between the piston and hydraulic cylinder.In this study, PTO stiffness and damping are 100,000 N/m and 320,000 N/(m/s), respectively.

Load boundary condition
Fig. 10 presents the load boundary condition used in this study.As can be seen from Fig. 10, the wave loads, which are obtained from the analytical model presented in Section 3.1, are applied to the arm tip.The direction of the load can be either vertically upward (positive load) or vertically downward (negative load).

Results and discussions
WECs are subject to significant cyclic loads induced by the wave, and therefore their design is generally dominated by fatigue.In this study, the MBD system is simulated under the fatigue DLC (design load case).In this case, the mean values of significant wave height H s and wave period T m at a specific spot in China's Bohai area from monitoring data provided by NOTC are used.The sea condition is presented in Table 4.
The wave loads under the DLC are calculated from the analytical model (see Section 3.1) and are presented in Fig. 11.From Fig. 11 we can see that 1) the wave loads vary periodically with the time; 2) the period of loading cycles for the DLC is 6.2s, which corresponds to the associated wave period.
Under the DLC, the results from the FMBD model presented in Section 4 are validated against the results from the analytical model.After the validation, the FMBD model is further used to calculate the stress distributions, fatigue life, deformations, modal frequencies and modal shapes of the structure, and acts as a reference point towards scaling and optimisation of the WEC.

Validation
The floater arm tip displacement and velocity calculated from the FMBD model are compared with the results obtained from the analytical model presented in Section 3.2.Fig. 12 depicts the comparison results.
From Fig. 12 we can see that 1) both vertical displacement and velocity vary periodically due to the periodic wave loads; 2) the results from the FMBD model show reasonable agreement with the results from the analytical model, both in terms of trend and magnitude, with a relative difference of 10.1% at the maximum value of the floater arm tip vertical displacement observed at the time of 6.2s and a maximum relative difference (5.6%) observed at the time of 2.2s for vertical velocity.This confirms the agreement between the FMBD and analytical models.
It is worth mentioning that the maximum relative difference of the vertical displacement between the FMBD and analytical models in this case is larger than that of the vertical velocity.This is probably due to the fact that the displacement is not only dependent on the velocity but also affected by other factors, e.g. the acceleration.Therefore, a small difference in velocity and acceleration may result in a relatively larger difference in displacement.

Applications
After the validation, the FMBD model is further applied to calculate the stress distributions, fatigue life, deformations, modal frequencies and modal shapes of the structure.

Stress distributions
Experimental fatigue testing data are mostly uniaxial, whereas FEA results are generally multiaxial.Von Mises stress is generally used to take account of the multiaxial stress state and to compare against the experimental uniaxial stress value.However, the stress range calculated from a strictly positive history of the von Mises stress does not include any potentially negative part of the stress cycle, limiting its application to fatigue assessment to some extent.
The signed von Mises stress s svm is thus used to mediate this shortcoming by applying the sign of the largest absolute principal stress s p to the von Mises stress s vm , i.e. [21]: The signed von Mises stress has been widely used for fatigue assessment [22e24], and therefore is chosen in this study for fatigue assessment.The time history of maximum signed von Mises stress of the floater arm, piston, hydraulic cylinder and bracket under the DLC is presented in Fig. 13.
As can be seen from Fig. 13, the maximum signed von Mises stress of all parts show periodic variations.The stress range level on the floater arm, bracket, hydraulic cylinder and piston show decreasing trend, with the highest stress range observed at the floater arm and lowest stress range observed at the piston.The peak positive von Mises stresses of all parts during the simulation are all observed at the times of 5.6s, 11.8s and 18s, which correspond to the peak negative wave load (see Fig. 11).The von Mises stress distributions of all parts at the time of 5.6s are presented in Fig. 14.As can be seen from Fig. 14 25.0 MPa and 32.2 MPa, respectively.The locations of the maximum von Mises stress, i.e. the hot points, are also presented in Fig. 14.

Fatigue life
A widely used method for fatigue analysis is the S À N curve method, which is based on fatigue testing data (i.e.stress range and associated number of cycles to failure).According to the S À N curve method, the number of loading cycles to fatigue failure, N, is given by: log where A is the intercept; m is the slope of the S À N curve in the loglog plot, DS is the stress range.
In this study, the values of both intercept A and slope m in Eq. ( 54) are taken from DNV standard [25] and are listed in Table 5.
The number of loading cycles, N, in Eq. ( 54) can be converted to fatigue life in year, N Year , by using the following equation: where T is the period of a loading cycle.
The stress range obtained from the FMBD model is then used in the S À N curve method to calculate the fatigue life.The maximum fatigue life is considered to be 20 years, which is identical to the design life.The fatigue assessment results of each part are presented in Fig. 15.
As can be seen from Fig. 15, the shortest fatigue life (2 years) is observed in the floater arm, and the longest fatigue life (20 years) is observed in the piston.Additionally, the fatigue lives of floater arm, hydraulic cylinder and bracket are lower than the target design life of 20 years.This indicates that the WEC prototype used in this study is prone to experience fatigue failure and needs to be improved in order to meet the target fatigue life.

Deformations
The time history of maximum vertical deformations of the floater arm, piston, hydraulic cylinder and bracket is presented in Fig. 16.It should be noted that the vertical deformations presented in Fig. 16 comprise both global displacements and local deformations.As can be seen from Fig. 16, the maximum deformations of all parts show periodic variations.The deformation levels of floater arm, piston, hydraulic cylinder and bracket show a decreasing trend, with highest deformation observed at the floater arm and lowest deformation observed at the bracket.The peak deformations of all parts during the simulation are all observed at the time of 3.6s, 10s and 16.2s.The deformations of each part at the time of 3.6s are presented in Fig. 17.In order to facilitate the illustration, the initial structure (i.e. the structure at time 0s) is also presented in Fig. 17.
As can be seen from Fig. 17, the maximum deformations on floater arm, piston, hydraulic cylinder and bracket are about 0.17 m, 0.068 m, 0.018 m and 7.87e-6m, respectively.The locations of maximum deformation are also presented in Fig. 17.As expected, the deformations of the bracket are negligible, due to the fact that the bottom surface of the bracket is fixed.

Modal frequencies and modal shapes
Modal analysis is used to obtain the dynamic properties of the MBD system, such as modal frequencies and modal shapes.In this case, the MBD system is considered to have free vibration (no loads on the structure).The first six modal shapes and frequencies obtained from the FMBD model are depicted in Fig. 18 and listed in Table 6, respectively.
As can be seen from Fig. 18, higher values of modal frequencies occur at higher modes, and different modal shapes associated with different modal frequencies show the pattern of vibration.Under the 1st, 2nd and 4th modes, the structure vibrates in the X-Y plane where the floater heaves.Under the 3rd and 5th modes, the structure vibrates out of the X-Y plane.The 6th mode is the torsional mode, in which twisted vibration occurs.
Wave-induced vibration is generally the main reason for the structural vibration of WECs, and therefore it should be considered in the design of WECs.In this case, the time period T m of the wave is 6.2s (see Table 3), and the associated wave frequency is 0.16 Hz.As can be seen from Table 6, the frequency of the lowest mode is 1.05 Hz, which is much higher than the wave frequency.This  indicates the structure is not likely to experience wave-excited vibrations.

Conclusions
In this study, an FMBD (flexible multibody dynamics) model for point-absorber WECs (wave energy converters) has been developed.The FMBD model is applied to the structural modelling of the NOTC (National Ocean Technology Centre) 10 kW multiple-pointabsorber WEC.The floater arm tip displacement and velocity obtained from the FMBD model are validated against the values from an analytical model, which is also developed in this work.After the validation, the FMBD model is further used to examine the stress distributions, fatigue life, deformations, modal frequencies and modal shapes of the structure.The following conclusions can be drawn from the present study: 1) Both the floater arm tip displacement and velocity obtained from the FMBD model show reasonable agreement with those from the analytical model, with a relative difference of 10.1% at the maximum value of the floater arm tip displacement, which confirms agreement of the FMBD and analytical models.
2) The frequency of the lowest mode is 1.05 Hz, which is much higher than the wave frequency.This indicates the structure is unlikely to experience wave-excited vibrations.
3) The deformation level of the floater arm, piston, hydraulic cylinder and bracket show a decreasing trend, with the highest  deformation observed at the floater arm and lowest deformation observed at the bracket.4) The stress range level on the floater arm, bracket, hydraulic cylinder and piston show a decreasing trend, with the highest stress range observed at the floater arm and lowest stress range observed at the piston.
5) Apart from the piston, the fatigue life of other parts (i.e.floater arm, hydraulic cylinder and bracket) of the WEC used in this study are lower than the target design life of 20 years, with the shortest fatigue life of 2 years observed in the floater arm, indicating that the WEC prototype used in this study is prone to experience fatigue failure and needs to be improved in order to meet the target fatigue life.
The FMBD model developed in this work is proven to be capable of the accurate modelling of point-absorber WECs, thus this approach can provide valuable information for designers to further optimise the structure and assess the reliability of structures of similar class.

Fig. 14 .
Fig. 14.Von Mises stress distributions at time of 5.6s: a floater arm, b piston, c hydraulic cylinder, d bracket.

Fig. 15 .
Fig. 15.Fatigue assessment results: a calculated fatigue stress range, b calculated fatigue life.

Table 2
Mesh independence study.

Table 4
Sea condition.