Modeling and Simulation of a Wave Energy Converter: Multibody System Coupled to Fluid-Film Lubrication Model and Thermal Analysis ^†

Abstract

Sea wave energy is being increasingly regarded as one of the most promising sources of renewable energy. This paper deals with the modeling and simulation of an onshore wave energy converter system designed by UMBRA GROUP SpA. Several topics are addressed. Starting from the multibody modeling strategy, this paper delves more deeply into the mechanical efficiency evaluation of the ball-screw in the elastohydrodynamic lubrication regime, the core of the energy conversion process, as well as the thermal characterization of the power take-off module, based on the lumped-parameter and finite element method models. High values of ball-screw indirect efficiency have been observed, ranging from 73% to 97%; these results appear even more encouraging when compared to the performance of alternative energy-consuming technologies. Thermal analysis, on the other hand, provided a maximum temperature increase of 40 °C, allowing for the aversion of any structural collapse and the realistic identification of the lubrication regime, which turned out to be mostly mixed. Finally, an inverse multibody dynamic analysis is performed, and the most interesting simulation results are collected to prove the effectiveness of the proposed approach.

1. Introduction

In the past few decades, the idea of exploiting the huge and largely untapped potential energy of ocean waves (typically 2–3 kW/m² compared to 0.4–0.6 kW/m² of wind and 0.1–0.2 kW/m² of solar radiation at the earth’s surface [1]) has drawn significant attention in order to meet renewable energy targets and environmental awareness. Great effort has been devoted to the development of several types of devices, commonly called wave energy converters (WECs), able to convert the immense power the ocean waves into electricity. According to the large energy density when compared with other energy systems, WECs can be a key player for increasing the sustainability of the energy sector [2,3,4,5,6,7].

Depending on their working principle, WECs are generally categorized into oscillating water columns, overtopping systems, bottom-hinged systems, and floating-point absorbers, which are thought to be the most cost-efficient technology to take advantage of energy density and the consistency of ocean waves. To date, companies and academic research groups around the world have modeled and tested many working concepts in wave tanks, but only a few prototypes have moved from the R&D stage to sea testing [8,9,10,11].

This paper deals with modeling and simulation of a point-absorbing type WEC, and the whole study is carried out in collaboration with UMBRA GROUP (UG), a leading global manufacturer of recirculating ball-screws in the aeronautics sector, which has been strengthening its R&D background in the energy market over the past ten years. As shown in Figure 1, UG developed a very promising concept for onshore WEC, consisting of a floating buoy connected to the support frame by a couple of holding arms. The buoy, acting as an energy absorber, is lifted and dropped by sea waves, oscillating in a single mode (i.e., it exhibits 1 DOF). Periodic forces exerted by ocean waves are transferred to the power take-off module (PTO), located between the floating buoy and the support frame. Within the PTO, a massive recirculating ball-screw directly converts the reciprocating movement of the nut into the high-speed rotary motion of the screw shaft, which can be easily fed to the electro-mechanical reciprocating generator (EMRG) for electricity generation.

Such a device is designed for shoreline installation, close to the utility network and it is expected to require less intensive maintenance compared to offshore WECs [12]. It is worth noting that, in many areas of the world, the wind blows with enough consistency and speed to provide continuous waves along the coast. However, as ocean waves are attenuated in their travel through the shallow water, the onshore location leads to a lower energy content available for regenerative purposes.

The aim of this work is to investigate the dynamic behavior of the WEC in common sea states by developing a multibody model supported by the analytical estimation of the ball-screw efficiency. Moreover, a standalone simulation of transient heat conduction within the PTO is carried out by relying on increasingly accurate models. Finally, the most interesting simulation results are collected, shedding light upon the effectiveness of the proposed models.

2. Multibody Model, Ball-Screw Efficiency, and Thermal Models

Primarily, it is wise to offer a focused, yet coherent, overview of the WEC multibody modeling strategy and the most interesting forces acting on the structure. As a first attempt, a simplified structure able to provide an overall idea of WEC functionality was developed in a multibody software environment, wherein all components were paired together by defining a local coordinate system coord and a set of constraints for each of them; the result is depicted in Figure 2.

At this stage, the rotor–stator interaction within the PTO is not yet considered, and its resistant action against bouncing motion of the buoy is performed by a cylinder-piston system, whose inner pressure p is supposed to be a function of rod position s and velocity v (1).

$$p=p\left(s\left(t\right),v\left(t\right)\right)$$

(1)

Subsequently, the first developed multibody model was moved to the MATLAB/Simulink environment, and the real structure of the WEC provided by UG was examined. In order to minimize the number of constraints to be imposed between the hundreds of disjointed parts (e.g., all nuts, bolts, screws, and levers), the 3D CAD provided by UG was preliminarily divided into five XML assemblies that behave like rigid bodies after being imported into MATLAB/Simulink: the onshore frame, the buoy and its connecting arms, the PTO shell and nut, the screw shaft and EMRG rotor, and the EMRG stator. Then, proper joints were included to impose kinematic constraints between the assemblies.

A multibody model is designed to receive either the vertical motion z(t) of the buoy CG, or the net force F_z(t) experienced by the floating body as inputs. For this purpose, both buoyancy force F_hs(t) and hydrodynamic force F_hd(t) have to be considered: the former is evaluated as the volumes of water V displaced by the completely or partially submerged buoy for each discrete angular position α, as further detailed in Figure 3, while the estimation of the latter requires more in-depth considerations concerning the velocity vector field of the flowing fluid around the buoy. Kinematic and dynamic inputs can be iteratively obtained by trying to achieve a match between the sea state energy E and the mechanical work W completed by the reciprocating nut and available for energy conversion, whereas “sea state” refers to the physics quantities {height H; period T} that fully define the sinusoidal sea wave. According to the linear wave theory proposed by G. B. Airy [13], the wave energy density E per unit of horizontal area is proportional to the wave height H squared (2), regardless of the depth of the water on whose surface the gravity waves propagate:

$$E=\frac{1}{16}\rho g{H}^2$$

(2)

where 𝜌 is the density of water, 𝑔 is the acceleration of gravity, and 𝐻 is known within physical oceanography as a significant wave height, but it can be assumed to be equal to the wave height, if more detailed wave characterization is not available.

The energy content transferred by the sea to the buoy can be obtained by multiplying the wave energy density E by an appropriate horizontal surface; common sense might suggest choosing as a surface reference the instantaneous intersection between the buoy and the free surface of the sea or, more conveniently, its average extension S_avg observed during usual WEC operating conditions. As a result, the sea wave lifts the buoy by an altitude h, starting from its static equilibrium position, causing the high-speed rotary motion of the screw shaft. The energy balance algorithm (EBA) represents the attempt to iteratively obtain the lifting h of the buoy and, ultimately, a kinematic input for the multibody model, starting only from sea state and the mandatory simplification of the hypotheses, as shown in Figure 4.

The axial force F_PTO resisting the relative motion between the translating nut and the rotating screw shaft is represented by a vector in Figure 5, and Equation (3) provides an estimation of it,

$$F_{PTO}=-K_{damp} v_{PTO}-F_{friction} sign\left(v_{PTO}\right)-K_{visc} v_{PTO}-K_{inert} a_{PTO}$$

(3)

where v_PTO and a_PTO represent the nut’s velocity and the acceleration of the translating nut, respectively; F_friction is the friction force, depending on which lubrication regime the recirculating ball-screw is expected to operate in; K_damp, K_visc, and K_inert are inertias and damping coefficients, as each subscript suggests. Since frictional forces have a large influence on this mechanism, the estimation of output power P_out available at the rotating screw shaft requires contact analysis. For this purpose, the multibody model described above is meant to be coupled to the elastohydrodynamic lubrication (EHL) model for the estimation of the indirect efficiency η of the recirculating ball-screw, in a simultaneous calculation scheme, according to Figure 6, with data transferred from the multibody module to the EHL module regarding current relative speeds and forces at the ball-screw mechanism level. As very well underlined in [2], the interdisciplinary research efforts between mechanical, electrical, and control engineers are required to achieve mature research in WEC systems. Thus, the involvement of a tribological coupling as the ball screw mechanism in the proposed WEC system requires an adequate investigation of the frictional response of the mechanism responsible for the motion conversion.

An accurate assessment of efficiency by relying exclusively upon the geometry of the contact surfaces and their finish, the helix angle of the thread, as well as the working conditions of the screw (load, speed, lubrication, preload, alignment) is challenging due to the complexity of the kinematic and dynamic analysis of the rolling elements. Indirect efficiency η is defined as the ratio between output power P_out available at the rotating screw shaft and input power P_in provided by the translating nut to achieve the conversion of motion against frictional resistance [12], as reported in Equation (4).

$$\eta =\frac{P_{out}}{P_{in}}=\frac{M_{scr}{\omega }_{scr}}{F_{nut}{v}_{nut}}$$

(4)

Since the frictional forces and the lubrication regime engaged during the operation have a large influence on the overall efficiency of this mechanism, the transmitted torque M_scr evaluation requires a contact analysis, whereas axial force F_nut and both velocities ω_scr and v_nut can be easily obtained from multibody simulations. The equilibrium condition in contact points A and B are depicted in the Figure 7, which shows the forces shared between the rolling element, the screw shaft, and the nut.

Let the following assumptions hold:

• pure rolling of the ball on the nut;
• constant slip angle γ_B at ball-screw interface.

Furthermore, in smooth and moderate operating conditions, low values of the screw’s angular acceleration pave the way for a huge simplification in the proposed dynamic analysis [14]: inertia forces can be neglected, and the quasi-steady-state assumption (QSSA) is also expected to be appropriate. Let Q_A = Q_B = Q denote the closing force, β_A = β_B = β the ball contact angle, α the lead angle, and f and ρ the coefficient of friction and the angle of friction, respectively. The force F_ball applied on each loaded ball and the transmitted torque M_scr are reported in Equations (5) and (6) [15].

$$F_{ball}=Q\left\{\cos \left(\rho \right)\left[\sin \left(\beta \right)\cos \left(\alpha \right)+f\cos \left(\beta \right)\cos \left(\alpha \right)\right]+sin\left(\rho \right)sin\left(\alpha \right)\right\}$$

(5)

$$\begin{gathered} M_{scr}=\left(r_m-r_b\cos \left(\beta \right)\right)Q\left\{cos\left(\rho \right)\left[\sin \left(\beta \right)\sin \left(\alpha \right)-f_B\left(\cos \left({\gamma }_B\right)\cos \left(\beta \right)\sin \left(\alpha \right) \\ -\sin \left({\gamma }_B\right)\cos \left(\alpha \right)\right)\right]-sin\left(\rho \right)cos\left(\alpha \right)\right\} \end{gathered}$$

(6)

Finally, the elastohydrodynamic lubrication occurring in nonconformal elliptical contacts is investigated [16,17] and the minimum film thickness h_min, as well as the coefficient of friction f and the angle of friction ρ, are predicted, employing the widely used Hamrock–Dowson formula [18], and the Houpert [19] and Balan [20] outcomes. It may be noted that the minimum film thickness h_min is often divided by the root mean square roughness σ of the contacting surfaces [21] in order to calculate the film parameter λ to predict the lubrication regime the machine is expected to operate under (7), as shown in Figure 8 [22].

$$\lambda =\frac{h_{min}}{\sigma }$$

(7)

Mechanical power losses taking place in the ball return system and supporting the bearings, as well as electrical power losses certainly confined inside the EMRG, cause the whole actuator to overheat. The decrease in oil viscosity due to higher temperature may result in a change in the lubrication mode from hydrodynamic to mixed film to boundary lubrication, which usually represents the critical condition that reduces the lifespan of the components. Moreover, the correct estimation of the EMRG operating temperature is helpful in order to prevent thermal damage and ensure its continuing operation. For this purpose, the lumped-parameter thermal model and the more advanced bidimensional thermal model based on the finite element method (FEM) have been developed for the offline estimation of temperature dynamics in the PTO module. Since heat production estimation deeply affects the time evolution of temperature, both mechanical and electrical power losses need to be taken into account. For a typical WEC operating cycle lasting 600 s, experimental data for instantaneous mechanical power losses are provided by UG and shown in Figure 9, whereas the electrical losses have been evaluated assuming an average generation efficiency of 87%, which is known to be the useful electric power output.

As first attempt, due to their relatively small sizes in the radial coordinate and their high thermal conductivities, the generator, the screw shaft, and the nut are supposed to behave as “lumps” whose internal temperatures remain essentially uniform in space. Under such an assumption, the temperature can be taken to be a function of time only [23], and the lumped parameter thermal model reported in the first order linear ODEs (8), (9), and (10) is carried out. Each component, initially at uniform temperature T₀, is suddenly subjected to heat generation $\dot{Q}$ and a convective boundary condition (i.e., convective heat flux condition into the surrounding environment). In this model, T_gen denotes the generator temperature, T_scr the screw temperature, T_nut the nut temperature, T_air,1 the temperature of the air around the generator, T_air,2 the temperature of the air near the nut and the screw, m the mass, c_p the specific heat capacity, ${\dot{Q}}_{el}$ the electrical power loss, ${\dot{Q}}_{mecc}$ the mechanical power loss, $a_1$ the convective thermal conductance referred to the generator-air interface, $a_2$ the convective thermal conductance referred to the screw–air interface, $a_3$ the convective thermal conductance referred to the nut–air interface, $b$ the contact thermal conductance for the screw–generator interface, and p stands for the partition coefficient of heat generation between the two components mentioned as a subscript (i.e., screw, nut, bearing) and assumed equal to 50% in both cases, while r is the fraction of mechanical losses taking place in the recirculating ball system supposed equal to 95%.

$${\left(m{c}_p\frac{dT}{dt}\right)}_{gen}=a_1\left(T_{air,1}-T_{gen}\right)+b\left(T_{scr}-T_{gen}\right)+{\dot{Q}}_{el}$$

(8)

$${\left(m{c}_p\frac{dT}{dt}\right)}_{scr}=b\left(T_{gen}-T_{scr}\right)+a_2\left(T_{air,2}-T_{scr}\right)+p_{scr-nut}r{\dot{Q}}_{mecc}+p_{scr-bea}\left(1-r\right){\dot{Q}}_{mecc}$$

(9)

$${\left(m{c}_p\frac{dT}{dt}\right)}_{nut}=a_3\left(T_{air,2}-T_{nut}\right)+\left(1-p_{scr-nut}\right)r{\dot{Q}}_{mecc}$$

(10)

The thermal parameters were obtained in a literature search [24,25] or, whenever the need arose, estimated using the parameter estimation tool in Simulink [26,27,28,29] by fitting the experimental profiles of thermal, mechanical, and electrical energies per cycle provided by the UMBRA Group.

The finite element method (FEM) is used to numerically solve the coupled conduction–convection problem, overcoming the strong limitations and the restricted range of applicability experienced in the lumped system formulation. The FEM thermal model developed in ANSYS required:

• laying out the physical structure;
• specifying the thermal properties of the materials and the internal heat sources;
• assigning the boundary and initial conditions.

In order to take advantage of the axial symmetry of both the geometry and the boundary conditions, the bidimensional structure depicted in Figure 10 is considered, and the adiabatic condition along the axis of symmetry y can be assumed.

3. Simulation Results

This section covers the most interesting findings from the simultaneous simulations performed by the Multibody and EHL models, as well as the comprehensive results obtained by means of the standalone thermal models.

3.1. Multibody and EHL Models

Inverse dynamic analysis is performed to better understand and improve the power takeoff (PTO) techniques and ultimately enhance the energy extraction from sea waves, bearing in mind that WEC dynamic behavior is closely related to the efficiency of energy conversion. In Figure 11, Figure 12 and Figure 13, a sinusoidal vertical motion z(t) lasting for the period T = 3.8 s of the common sea states, with an amplitude A = 1 m, is provided as kinematic input to the buoy CG, and three outputs are observed: axial force F_PTO acting on the translating nut, screw rotation speed ω_screw, and nut velocity v_nut.

The EHL model provided film parameters λ and coefficients of friction f, shown in Figure 14, and referred to each component mentioned above.

Finally, indirect efficiency η is reported below in Figure 15.

The ball-screw frictional behavior used as the input of the lubrication model for the calculation of the indirect efficiency η results in high values of the latter variable, ranging from 73% to 97%, among the systems showing efficiency in the highest ranges, as explained in [30].

Moreover, the lubrication model proves that the ball-screw mechanism operates largely under mixed-mode lubrication, in the right-hand range of the Stribeck curve, associated with a decrease in the coefficient of friction f as the relative velocity increases.

3.2. Thermal Models

The lumped-parameter model has been implemented in MATLAB/Simulink, and time-domain solutions have been achieved in the offline mode. The simulation scenario is:

$$T_0=20 °\mathrm{C}, {\beta }_{scr-nut}={\beta }_{scr-bea}=50\%, \gamma =95\%$$

(11)

The geometrical and thermal properties of the components have been provided by UG, while the thermal conductance was obtained from the literature. In Figure 16, numerical solutions (dashed line) and approximating exponential functions (solid line) on the type $T\left(t\right)=T_0+\mathsf{\Delta }T\left(1-e^{-t/\tau }\right)$ are depicted. The careful choice of simulation times equal to 100 min allows the temperature to reach the steady-state value in each component.

Electrical losses four times higher than frictional losses result in lower thermal stress on the nut, whereas the screw and the electro-mechanical reciprocating generator experience a more noticeable, but not critical, temperature increase: simulations show that the maximum excursion between the warmer and the colder control points is only a few tens of Celsius degrees. Finally, the nut exhibits a uniform temperature increase just above 40 °C.

The FEM thermal model developed in ANSYS provided the most accurate thermal profiles: the spatial distribution of the temperature within the most thermally stressed component, as well as its time evolution in some control points, are shown hereafter in Figure 17 and Figure 18. As result, the figures show the compatibility of the temperatures reached by the main components of the WEC with the material properties and their resistance to thermal stress.

4. Concluding Remarks

This multidisciplinary research focused on the development of strongly integrated models to simulate dynamic, tribological, and thermal responses of a wave energy converter and its energy conversion process. As matter of fact, when looking at the present literature, one finds that the researchers who are focusing on the mechanical aspects of the wave energy converter systems and the device responsible for the energy capturing are often working apart from considering the behavior of electrical systems and vice versa, including the relevant coupling between thermal field. Thus, an effort to deepen the main frictional conjunction responsible for the essential stage of the energy conversion in this WEC system has been proposed in the present paper.

Inverse dynamic analysis demonstrates that high rotational speeds can be easily obtained, exploiting shallow water wave energy. It also provides the internal force F_PTO, evenly distributed among the rolling elements of the ball-screw and used as the input of the elastohydrodynamic lubrication model for the online estimation of the indirect efficiency η: high values ranging from 73% to 97% have been observed. Moreover, the lubrication model proves that ball-screw mechanism operates largely under mixed-mode lubrication, in the right-hand range of the Stribeck curve, associated with a decrease in the coefficient of friction f as the relative velocity increases.

Two thermal models have been compared to deliver additional results to designers of wave energy converters based on the ball-screw mechanism. Electrical losses four times higher than the frictional losses result in lower thermal stress on the nut, whereas the screw and the electro-mechanical reciprocating generator experience higher temperatures, operating in safe conditions in all the simulated scenarios.

For future study, the improvement of the design of the mechanical element will aim at proving the stability of such systems for the sea environment, including the rheological properties of greases and lubricants exposed to sea water. The design of the system should also encompass the automatic control of auxiliaries to manage system operating conditions outside of those reported in this work.

Researchers in the field wave energy converters will also be focused on the optimization of the energy transfer from the mechanical device to the grid on the other side in order to increase the capture of energy in a broad range of sea states.