Mechanical design and modeling of a single-piston pump for the novel power take-off system of a wave energy converter

A multi-pump, multi-piston power take-off wave energy converter (MP2PTOWEC) has been proposed for use with a novel renewable energy harvester termed the Ocean Grazer. The MP2PTO WEC utilizes wave motion to pumpevia buoys connected to pistonseworking fluid within a closed circuit and store it as potential energy that can be converted to electricity via turbines. This paper introduces the mechanical design and model-based performance prediction of a single-piston pump that constitutes the basic building block for the MP2PTO WEC. Results provide preliminary validation of aqueous lubrication as a viable means of reducing friction and wear, suggesting that water-based hydraulic fluids can prohibit solid contact at the piston-cylinder interface while reducing volumetric leakage, and allowing for an estimation of the energy extraction efficiency for the mechanical pumping system. Pending more thorough and extended tribological investigations using the methodology introduced in this paper, findings suggest that the overall system efficiency will be dictated by the hydrodynamics of the buoys actuating the pumping system. © 2016 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
A number of near-and off-shore wave energy converters (WECs) have been proposed in recent years based on attenuator [1,2], point absorber [3], overtopping [4], and other design principles; for a comprehensive review of existing WEC technologies, the reader is referred to a recent report published by the Strategic Initiative for Ocean Energy [5]. In an effort to improve on the state-of-the-art, the University of Groningen has patented a novel semisubmersible renewable energy harvester, termed the Ocean Grazer, with a multi-pump, multi-piston power take-off (MP 2 PTO) WEC at its core. A single Ocean Grazer device, for which the MP 2 PTO WEC employing multiple multi-piston pump units will contribute about 80% of power generation (secondary technologies such as oscillating water column and wind turbine systems will contribute the rest), is projected to produce more than 200 GWh/ year and have a storage capacity of about 800 MWh [6,7].
The operating principle of the MP 2 PTO WEC, shown in Fig. 1(a), is to create pressure difference (hydraulic head, H) in the working fluid circulating between two reservoirs that can be transformed into electricity via a turbine (T). A modified point absorber design will be used such that floating buoys (B i ) will follow the motion of an incident wave and actuate linear hydraulic pumps (P i ) to move the working fluid column during the upstroke. In this manner, the working fluid can be pumped to the upper reservoir where it will be stored as lossless potential energy, and allowing for the decoupling of electricity generation from the variability of available wave energy over timescales of seconds (for individual waves) to hours and days.
Varying sea conditions determine an incident wave's characteristics such that each wave may differ significantly from those preceding or following it; ideally, a WEC should be able to extract energy from both small and large waves with a range of periods. Energy extraction is expected to diminish the energy content and height of a wave as it moves through a WEC. In the case of the Ocean Grazer, which will employ a grid of buoys e termed a floater blanket e to actuate the pumps, the first pump unit can potentially extract more energy than the second, and so on. To account for the inherent variability in wave energy content, multiple pistons will be activated within each pump in the Ocean Grazer, as shown in Fig. 1(b), to maximize energy extraction for waves ranging in height from 1 to 12 m and periods of 4 to 20 s [6,7]: controlling the coupling between any buoy (B i ) and a number of variable-size pistons (P i,j ) can optimize the load the buoy has to carry to achieve resonance during the upstroke and let an uncoupled piston sink due to its own weight during the downstroke. Preliminary work has demonstrated the successful potential use of variableload control for a multi-piston pump [8]; however, the focus of the present work is the tribological characterization of the behavior of a single-piston pump unit e that will constitute the basic building block for the multi-piston and, eventually, the MP 2 PTO WEC e in order to better understand the dynamic behavior of the piston based on physics-based formulations of the forces and pressures acting on it. While variable-load control is not discussed within the context of the present model, this work will be essential in the future design and implementation of improved control schemes [9]. Similarly, previously developed methodologies for a hydro-pumped storage system utilizing a number of pumps in parallel operation will be useful in the eventual transition from a single-to a multi-pump system [10,11].
The main tribological interfaces of interest in the MP 2 PTO WEC are between the piston and cylinder and at the seals that isolate the working fluid from sea water while allowing for the cable connections to transfer the buoy motion to the pistons. While the sealing problem is important because of sea water's highly corrosive and biofouling attributes, it is not addressed in the current work. Instead, it is assumed that perfect sealing is achieved so that the working fluid circulation is completely isolated from sea water, while buoy motion is transferred without frictional losses to the piston via a cable or rod of known stiffness. Future research on the tribology of the cable-seal interface will allow for the relaxation of these assumptions, especially within the context of flexible diamond-like carbon coatings that can be used in rubber seals [12]. Flexible coatings will also be relevant in secondary tribological interfaces between the piston and piston flaps as well as the oneway valves preventing working fluid backflow.
Scaling and maintenance issues require the design and use of a robust and abundant lubricant supply to the tribological interfaces. This suggests that a water-based hydraulic fluid optimized for the operating conditions of the piston-cylinder interface could potentially be used as a lubricant. Aqueous lubrication [13] was initially investigated within the context of environmental humidity at tribological interfaces [14], and subsequently for use with polymer composite coatings [15e17], and ceramic coatings with [18] and without texturing [19]. As a starting point, the present work assumes that a water-based hydraulic fluid (e.g. ISO Class HFC, a watery polymer solution) is used as the working fluid, while thermal effects are neglected at the piston-cylinder interface, where elastohydrodynamic (EHL) lubrication regimes are expected. Similarly to reciprocating engines [20], Sterling engines [21] and ring-less compressors [22], piston motion could result in boundary lubrication when the relative velocity at the interface becomes close to zero (at the top-and bottom-dead-centers), while EHL can be maintained otherwise. One of the questions to be answered with this work is whether EHL can be maintained throughout the piston stroke. Focusing on the piston-cylinder interface, the present model adopts existing methodologies to analyze the incompressible, inviscid and isothermal EHL problem. Future work will investigate the potential inclusion of piston rings, which may reduce working fluid backflow but increase friction, thereby complicating the lubrication issue.
A two-degree-of-freedom (2DOF) dynamical model comprising switching-state (upward versus downward motion), second order equations with an external forcing function was formulated in previous work [7,8]. This did not include friction and simplified the lubrication regime at the piston-cylinder interface by assuming hydrodynamic (Couette) flow. The current model builds on this dynamical model by adding EHL and solid contact and friction forces at the piston-cylinder interface based on solutions of the EHL problem. Further improvements to the model are planned in future work, including the formulation of realistic (irregular) wave displacement profiles serving as the external excitation, and accounting for the loss of wave energy content (and height) in the downstream direction due to energy extraction relevant for the MP 2 PTO WEC; the inclusion of additional degrees of freedom for the piston and the buoy, such as downstream translation depending on buoy hydrodynamics; and, the improved modeling of dynamical piston behavior including piston flap (valve) dynamics and the drag force acting on the piston during sinking within the fluid column as a function of piston design.
The overall efficiency of the Ocean Grazer will depend on frictional and hydrodynamic energy losses at the MP 2 PTO WEC. Therefore, minimizing friction and understanding the hydrodynamic behavior of grids of buoys are of paramount importance to maximize wave energy extraction as is reducing wear to ensure robust operation with little need for maintenance, especially in the presence of multiple tribological interfaces. Our results show that volumetric leakage at the piston-cylinder separation is critical in determining the overall pumping efficiency: this will decrease to below 80% for piston-cylinder separations larger than 200 mm, a finding that agrees with computational fluid dynamics simulations [23], pointing to target piston-cylinder separations around 100 mm when pure water is used as a lubricant. Such separations are comparable to those used in relevant applications such as ringless compressors, diesel and Sterling engines with typical values of 10 mm, 63 mm and 500 mm, respectively, and their optimization will be the focus of future work [20e22]. The present model is the first step in validating the pumping system's mechanical efficiency that has been measured in a small-scale prototype to be close to 99% and, hence, the potential of the Ocean Grazer to becoming a viable renewable energy harvester.

Piston excitation due to wave motion
The issue of wave hydrodynamics and their effects on floating structures is not trivial and has received significant attention in the literature [24,25], especially with reference to the control of point absorber WECs [26,27]. The scattering of an incident wave induces vertical (heaving) buoy motion, which results in waves radiating away from the oscillating buoy; in turn, momentum changes in the surrounding fluid give rise to net forces acting on the buoy. When isolating a single floating member for analysis as we do in this work, we choose to neglect the effects of neighboring buoys to the wave motion, energy content and resulting forces. In this section, we present the equation of motion e including the effects of buoyancy, wave excitation and drag forces e for a single prismatic buoy connected to a piston via a rod or cable of known stiffness under the excitation of a simple harmonic wave. The proposed model utilizes numerical methods and can therefore be used to predict the system response to more complex wave profiles, which can be characterized a priori as time histories of wave surface displacement z w (t). Furthermore, the dynamic adaptation of wave energy content with energy extraction, which will affect downstream wave height, is not relevant for a single buoy and is neglected in the present model; however, this will be included in the future modeling of the MP 2 PTO WEC that will utilize grids of buoys (a floater blanket).
Let an incident wave have a sinusoidal shape with its displacement about a zero mean, velocity and acceleration defined as where H w is the wave height and T w is the wave period [24]. At any time instant, a prismatic buoy of mass m b , height H b and crosssectional area A b will have displacement z b and will be submerged in the sea water by an amount D b defined as Consequently, the buoyancy force acting on it will be where r sw is the density of sea water and g is the gravitational constant. The buoyancy force will be constant for a fully-submerged buoy, zero if the buoy comes out of the sea water, and dependent on the magnitude of D b otherwise. When the buoy moves in stationary fluid, its acceleration induces the fluid in its immediate neighborhood to accelerate and this, in turn, induces an added mass effect onto the member [24]. This added mass can be described in terms of an added mass coefficient C a [25], and varies with the amount of submersion. Accounting for the added mass, the equation of motion of the buoy connected to a piston with relative displacement z and velocity _ z takes the following form: where K(z b Àz) ≡ F K and Cð _ z b À _ zÞ≡F C are the spring and damping forces representing the connecting rod (or cable) between the buoy and piston, while F d and F e represent the drag and excitation forces, respectively. The value of the damping ratio between the buoy and piston can be used to calculate the damping coefficient , with masses m 1 and m 2 defined in a later section (x2.4). The connecting rod/cable has area A r , density r r and elastic modulus E r so that its stiffness can be expressed as K ¼ E r A r /L r while its mass is m r ¼ r r A r L r . The weight of the rod will be accounted for in the piston equation of motion.
The drag force resists the motion of the buoy and is a function of the buoy velocity and drag coefficient C d : Published data can be used to extract a steady drag coefficient C ds for steady flow as a function of buoy geometry [25]. The drag coefficient C d in Equation (6) is the product of this steady drag coefficient C ds with an amplification factor quantifying the relative effect of inertial versus drag forces. In turn, the amplification factor is a function of the Keulegan-Carpenter number that depends on the maximum buoy velocity, the wave period and the width of the prismatic buoy. The KCnumber for the proposed system will be in the range of 0 < KC < 10, meaning that the drag coefficient can be several times larger than the one obtained from steady flow data [25]; however, we make the simplifying assumption that C d z C ds and plan to investigate this issue via experimental and simulation studies in future work.
The excitation force of Equation (5) includes pressure, inertial and damping contributions: where the buoy submersion has been defined in Equation (2) and the wave number is k w ¼ 2p/l w for wave length l w [28]. For the purposes of the present work, and in the absence of relevant data, we assume that B ¼ 0 in the current simulations. This is an important assumption since this damping coefficient will affect the dissipation of wave energy and, hence, the energy available for extraction. Preliminary investigations suggest that the magnitude of the wave damping coefficient is~35 Â 10 4 Ns/m for the prismatic buoys proposed in the mechanical design of the Ocean Grazer. Simulations adopting this value for the wave damping coefficient yield very comparable results to those reported herein for B ¼ 0; hence, we chose to not include them in this report pending more accurate hydrodynamic analyses. Similarly to the drag coefficient issue, we plan to perform experimental and numerical work to characterize the damping response of the buoys to be used in the full system for various buoy designs. This information will be used in future iterations of the model.
It should be noted that the conditions where the buoy comes out of the sea water completely or is fully submerged do not occur in our simulation results; instead, the buoyancy and excitation forces continuously scale with the instantaneous submersion calculated from the imposed wave displacement and the buoy heave motion as defined in Equation (2).

Pumping force and working fluid conservation
Following the buoy displacement, a cylindrical piston of height H p , radius R p and mass m p moves within a cylinder of effective length L c and cross-sectional area A c to pump the working fluid of density r from a lower to an upper reservoir. The effective cylinder length is measured relative to the bottoms of the two reservoirs, as shown in Fig. 2. The piston is hollow so that the top piston surface separates the two control volumes of working fluid with heights L 21 and L 43 . The top piston surface comprises two flaps (e.g. of semicircular shape) that can open and close passively with the flow. Check valves control the flow from the cylinder to each reservoir: during the upstroke, the piston flaps close and the check valves open, allowing working fluid to enter the cylinder from the lower reservoir and discharge into the upper; conversely, during the downstroke, both check valves are closed and the piston flaps open, allowing the piston to sink into the stationary fluid column. Flow within the cylinder of cross-sectional area A c is channeled from and to the upper and lower reservoirs with cross-sectional areas A U and A L , respectively. The initial values of the hydraulic heads (fluid levels) within the reservoirs are denoted as L U and L L .
As an initial condition, let the piston center of mass (COM) be located at an initial height L 0 measured from the bottom of the lower reservoir: this location will correspond to the center of the piston stroke relative to the cylinder base and should therefore be larger than half of the wave amplitude. Since the top piston surface thickness is very small relative to the fluid column, it can be neglected in the definition of L 21 and L 43 . Then, the dynamically varying distance between points 4 and 3 ( Fig. 2) is where L COM is the offset between the COM and the top surface of the piston. Similarly, the distance between points 2 and 1 is The magnitude of the force required to pump the fluid column of the upstroke, hereafter termed the pumping force, is where points 2 and 3 correspond to the top and bottom surfaces of the piston flaps and the cylinder cross-sectional area, depending on the piston radius R p and the piston-cylinder separation s, is The flow rate of the piston pump is a function of piston velocity Capacitance is derived from working fluid volume conservation and relates the pressure changes of the upper and lower reservoirs to hydraulic head fluctuations in each reservoir. Since flow rate is defined as the change in working fluid volume (area Â height), the change in hydraulic head can be expressed as _ H ¼ Q =A, and the change in pressure at points 1 and 4 becomes The change in the hydraulic head of the lower reservoir is negative as the working fluid flows from the lower to the upper reservoir during pumping. Integrating equation set (13) will yield the instantaneous values of the pressures at the upper and lower reservoir p 1 and p 4 , respectively, which serve as state variables in the model.
Continuity in the working fluid system requires that the following equation holds: (14) so that the pressure gradient between points 2 and 3 can be expressed as The pressure gradients between points 2-1 and 4-3 can be calculated as functions of the inertance, which is a measure of the pressure difference required to cause a change in flow rate with time, and the dynamically varying head (neglecting losses due to wall friction, etc). Specifically, the pressure gradients are (16) with the change in flow rate calculated by differentiating Equation (12) with respect to time the dynamically varying hydraulic head in the upper reservoir calculated as a function of the pressure at point 1, i.e. L U ¼ p 1 /rg, and the term r _ z 2 representing the change in momentum needed to accelerate the working fluid that enters the cylinder from the lower reservoir to attain the velocity of the moving fluid column. This formulation assumes that a check valve separates the working fluid volumes in the upper reservoir (V U ¼ A U L U ) and the cylinder (V c ¼ A c L c ) during the piston downstroke; however, this valve will open passively in the upstroke to accelerate the working fluid column whose total height becomes L 21 þL U . Since Combining the above formulae into the expression for the pumping force yields where the second term includes the inertia and weight of the fluid column of nominal mass m f,nom . It should be noted that cavitation is expected to occur for tall fluid columns. Placing the center of the piston stroke close to the bottom of the cylinder should reduce negative pressures during the upstroke: this is realized by appropriately selecting the value of L 0 to be small relative to the total cylinder length L c but still larger than H w /2. In addition, pressurization of the reservoirs, achieved by applying additional pressures p L and p H as parts of the initial conditions of the system, can also be used to avoid cavitation. Presently, we assume that p L ¼ p H ¼ 0.
The pumping force becomes zero during the downstroke for v z < 0. This may introduce discontinuities to the system since the inertial term m f ;nom ð€ z þ gÞ can be very large in the upstroke and results in an impact force (shock) acting on the piston during switching between the up-and downstrokes if this occurs instantaneously. In the present model, exponential growth and decay terms are introduced to the calculation of the fluid column mass m f to compensate for this behavior: where G is a growth/decay factor, t up denotes the starting time of the upstroke while m up,final and t up,final correspond to the fluid column mass and time instance at the last iteration of the upstroke, respectively. In addition, the value of m f is bounded by 0 and the nominal value m f,nom defined as and derived in Equation (18). A more nuanced description could introduce a flow-dependent function for F p that goes to zero more gradually as the piston flaps open incrementally in the downstroke based on predictions of their dynamic response. This, along with other issues such as the resistance of fluid flow within the cylinder due to wall friction, will be addressed in future work with the aid of computational fluid dynamics simulations. The authors have performed similar work that will be relevant to the unsteady flow in the section below the piston, in order to compute the exact velocity and pressure field and the minimum developed pressures [29,30].

Elastohydrodynamic lubrication at the piston-cylinder interface
Fig . 3 shows the side (a) and top views (b) of the piston-cylinder interface with the parameters relevant to the formulation of the elastohydrodynamic lubrication (EHL) problem [31]. For the purposes of this study, the tilting of the piston is only accounted for in the solution of the EHL problem and the tilting angle is assumed to be very small otherwise. The connecting cable or rod is attached at the piston COM, located at a distance L COM below the top surface of the piston of radius R p and height H p . The piston is contained within a cylinder of radius R p þs, where s is the separation between the piston and cylinder surfaces. Eccentricity e is defined as the radial offset between the COMs of the piston and cylinder. Assuming that the connecting cable/rod of length L r can be represented as a straight line, an angle 4 is formed between the tilted piston and the z-axis so that the eccentricity can be calculated as e ¼ L r sin 4zL r 4 (21) for very small values of 4. A cylindrical coordinate system (r,q,z) with its origin at the piston COM is adopted for the solution of the EHL problem, as shown in Fig. 3. Taking advantage of axisymmetry, the corresponding Cartesian coordinate triad (x,y,z) can be rotated about the z-axis so that the eccentricity is always defined along the x-direction. Therefore, tilting of the connecting cable during the upstroke creates a wedge in the fluid film with the minimum film thickness occurring at the bottom and top of the piston during the up-and downstroke, respectively. The wedge effect results in the localized buildup of fluid film pressure that may induce elastic deformation in the surrounding solid surfaces [31]. For simplicity, and since piston tilting is accounted for in the EHL problem alone, EHL forces are applied at the piston COM so that moments about the y-axis are neglected in the formulation of the dynamical model. Indeed, while the maxima of the hydrodynamic forces are expected to occur around the location of the minimum film thickness, the film pressure is distributed over the piston surface so that the abovementioned assumption should hold for very small tilting angles.
Given that the origin is taken at the piston COM, the center of the cylinder can be thought of as being offset by an amount e in the negative x-direction for a positive (counterclockwise) tilting angle 4, defined in the third quadrant as shown in Fig. 3. Tilting also results in the shifting of the piston COM by an amount dz ¼ L r (1Àcos4); this can be neglected for a very small tilting angle since cos4 z 1. The detailed formulation of the EHL problem, including the lubricating film thickness, can be found in Appendix 1.
The two-dimensional Reynolds equation in cylindrical coordinates, based on the incompressible and inviscid Navier-Stokes equations and the continuity condition [22,31] can be solved to yield the value of the film pressure p(q,z). A full film is assumed to exist at all times within the radial clearance over the range of q and the height of the piston. Fluid flow is assumed to be incompressible, inviscid with dynamic viscosity m, and laminar since the piston-cylinder separation is significantly smaller than the piston dimensions, while pressure variation in the radial direction is neglected for the same reason. Piston velocity is calculated from the dynamical model (refer to Fig. 2) as v z ¼ _ z. The following boundary conditions are assumed to hold:  (23) and are used to solve the Reynolds equation. In addition to the Neumann and Dirichlet boundary conditions specified in equation set (23), the Half-Sommerfeld boundary condition is also implemented whereby pressures smaller than the vapor pressure of the working fluid (absolute vapor pressure p v ¼ 2.2 kPa for an ISO Class HFC water-based hydraulic fluid at 20 C) are set equal to the vapor pressure [31]. Nevertheless, cavitation does not appear in the simulation results due to the selected location of the piston stroke center close to the bottom of the cylinder.
The values of the pressures can be calculated from equation set (16) as functions of the state variables p 1 and p 4 during the upstroke: under the assumption that the piston height is much smaller than the cylinder height and including the effect of leakage through the piston-cylinder interface (note: L U ¼ p 1 /rg). Instantaneous leakage is given by [22] Q leak ¼ and the corresponding leakage velocity is During the downstroke, which occurs when the criterion v z < 0 is satisfied, the pressures at the top and bottom of the piston flaps equalize. Flow into and out of the cylinder is inhibited by check valves; this is imposed by maintaining p 1 and p 4 constant, Q ¼ 0 and setting the fluid column acceleration to zero. Hence, pressure at the instantaneous piston location can be calculated as: EHL at the piston-cylinder interface will result in two fluid forces acting on the piston: a normal (bearing) force pushing the piston towards the center of the cylinder and a shear force resisting the vertical piston motion. These can be calculated using the following equations [22,32]: A penalty method is used to account for solid contact whereby an elastic (Hertzian) restoring force is activated whenever the fluid film thickness becomes negative [33]: where x L ¼ xþ(H p ÀL COM )sin4 is the largest lateral displacement measured at the bottom of the tilted piston. In potential instances of solid contact, the film pressure and thickness distributions are maintained constant and equal to the last step where EHL existed. Again, the reader is referred to Appendix 1 for descriptions of the variables used in the EHL model.

Equations of motion for the full system
The buoy of mass m b tracking the wave displacement z w displays heaving motion z b and actuates a piston with displacements x and z in the lateral and vertical directions respectively, as shown in Fig. 4 (also refer to Figs. 2 and 3). As stated previously, tilting of the piston is accounted for only in the solution of the EHL problem and is assumed to be negligible in the dynamics of the system (i.e. moments about the piston COM are neglected); however, the force acting to restore the piston COM to the center of the cylinder is accounted for. Displacements (DOFs) are measured relative to the COMs of the buoy and piston, respectively. Free body diagrams for the up-and downstroke are shown in Fig. 4: the pumping force is acting only during the upstroke and becomes zero otherwise.
In the case of the buoy (not shown in Fig. 4), the force balance in the vertical direction has been derived in a previous section. If we define an equivalent buoy mass Equation (5) becomes Referring to the free body diagram of the piston in Fig. 4(a), the force balance in the vertical direction has the following form during the upstroke: Introducing the pumping force of Equation (18) into the above expression yields: Let us define the equivalent mass m 2 in the following manner: Fig. 4. Free body diagrams of the piston in the upstroke (a) and downstroke (b).
in the upstroke in the downstroke (35) with the fluid column mass having been introduced in Equation (19). The equation of mass of the piston then becomes with the fluid mass m f becoming zero and the pressures p 1 and p 4 being kept constant during the downstroke. The force F f gives way to the solid friction force m st F r in the case of solid contact.
The force balance equation in the lateral direction is where the normal force F n acts during EHL and becomes the restoring force F r in the eventuality of solid contact. The force acting to restore the piston COM to the cylinder center for a small (positive) tilting angle 4 is derived from the equations of motion of a pendulum. The force balance in the lateral direction can therefore be written as follows: In addition to Equations (31), (34) and (39), equation set (13) must be solved simultaneously to yield the instantaneous values of the pressures p 1 and p 4 . In matrix form, the system equations follow the general form where the state vector is The matrix A is and the forcing vector takes the following form: where m st is the friction coefficient for steel-on-steel contact (assumed to be constant). Different materials and coatings are expected to be used in the actual system; hence, the friction coefficient will be different, while more accurate models of contact and friction can also be used without needing to define a friction coefficient a priori [34]. During the downstroke, check valves are activated to prevent backflow and the piston is allowed to sink passively into the working fluid column with open flaps. This is accounted for in the model by maintaining the pressures p 1 and p 4 constant and equal to the values attained at the last iteration of the upstroke. The ODEs of Equation (40) are solved numerically using the Runge-Kutta-Fehlberg algorithm with adaptable step size [35]. The initial conditions of the state variables for the full system are based on the initial wave surface position z w,0 with initial lower reservoir pressurizations p H,L and initial equivalent mass values of The initial displacements of the buoy and piston are calculated based on their equilibrium positions and accounting for the buoyancy force. It should be noted that the vertical displacements of the buoy and piston are expressed relative to a zero equilibrium position but will be separated by the connecting rod length L r in the physical system.

System parameters and assumptions
Simulations were performed for a single-piston pump with the resulting time histories of displacements, forces, flow-rates and pressures discussed in the following sections. Simulation parameters are summarized in Table 1. Table 2 summarizes material and other constants, while Table 3 lists the parameters used in the generation of a harmonic wave. The working fluid is assumed to be a water-based hydraulic fluid (ISO Class HFC) at 20 C and any temperature or pressure effects are neglected. Nevertheless, we can reasonably assume that the temperature would remain close to the nominal value (not necessarily 20 C) due to the flooded conditions of operation. The average density of HFCs, adopted from the specifications given for HFC by Bosch Rexroth (document RE 90233, January 2015), is assumed to be 1.08 kg/cm 3 , while the kinematic viscosity based on the ISO VG 68 specification is n ¼ 68 cSt; therefore, the absolute (dynamic) viscosity can be calculated to be m ¼ rn ¼ 0.0734 Pas. Such waterbased hydraulic fluids (containing up to about 98% vol% water) have been shown to behave as Newtonian fluids [36], further suggesting that temperature or pressure effects can indeed be neglected. In comparison, water at 20 C has the following properties: r w ¼ 998 kg/m 3 , m w ¼ 1.002 mPas, and p v,w ¼ 2.3 kPa, i.e. its viscosity is roughly 70 times smaller than that of HFC. Water was used in additional simulations and juxtaposed with the HFC results, while a piston-cylinder separation s ¼ 100 mm was selected to reduce volumetric leakage. Simulations performed with a commercial finite element solver (COMSOL Multiphysics software) have shown that, with water at 20 C as a lubricant, the pumping efficiency decreases to below 90% when the piston-cylinder separation exceeds 250 mm, and becomes smaller than 20% for separations of 1 mm [23]; this is a result of the high volumetric leakage under high pressure differences that increases with increasing separation. The higher viscosity of HFC allows for larger piston-cylinder separations able to maintain a high pumping efficiency while simultaneously reducing the probability of solid contact, as will be discussed below.
The structural parameters of the buoy, connecting rod/cable, piston and reservoirs are listed in Table 4. Damping between the buoy and piston is characterized by the damping ratio that is arbitrarily assumed to have a value of 5%, while no initial pressurization was applied to the reservoirs since the pressure of the working fluid in the system never decreases below the relative vapor pressure. The initial transients corresponding to the first period (0e10 s) are the result of the initial conditions applied to the state variables during the numerical solution; these have been removed from the plots for clarity and, as a result, the time histories now span the range from 10 to 50 s. Even in the absence of wave damping (B ¼ 0), the buoy is shown to smoothly track the wave motion after some initial transient. Nevertheless, an accurate characterization of this damping coefficient, possibly including linear and quadratic terms [25], will be necessary for the implementation of more complex wave profiles. The net vertical force acting on the buoy is given by

Time histories of displacements, forces, flow-rates and pressures
where the buoyancy, drag and excitation forces are defined in Equations (3), (6) and (7) and plotted in Fig. 5(b). The net force acting on the buoy will be transmitted to the piston via the connecting rod/cable whose design will in turn affect the dynamic response of the piston. Piston velocity serves as the criterion for determining upward versus downward motion as shown in Fig. 6(b): v z > 0 during the upstroke and v z < 0 during the downstroke. As would be expected, this corresponds to changes in vertical piston motion as shown in Fig. 6(a) and denoted by the overlay plots of a 'switch' variable that becomes one during the upstroke and zero during the downstroke.
Switching between the up-and downstrokes gives rise to some higher-frequency oscillations observed in the velocity and shear force time histories plotted in Fig. 6(bed), in addition to the plots of the 'switch' variable, due to the effect of the fluid column mass that goes from zero to several thousand kilograms or vice versa, resulting in an effective impact force acting on the piston at the instances of switching. Indeed, the frequency of these oscillations (~5.5 Hz) correlates well with the fundamental frequency of the connecting rod/cable f ¼ (1/2p)K/m 2 during the downstroke (~4.8 Hz).
Water at 20 C can initially maintain EHL at the piston-cylinder interface without solid contact, which begins with the second cycle, as shown in the insets of the shear and fluid bearing (normal) forces in Figs. 6(d) and 7(d); bouncing is observed inside the cylinder walls once solid contact is initiated. Since solid contact is expected to result in wear of the piston and cylinder surfaces, working fluid viscosity should be increased at the piston-cylinder interface in order to minimize solid contact. This is clearly shown in the main results of Figs. 6(aec) and 7(aec) pertaining to the HFC simulations: the pressure buildup in the lubricant is sufficient to separate the solid surfaces without solid contact, while the increased lubricant viscosity allows for a larger piston-cylinder separation (s ¼ 400 mm for HFC compared to 100 mm required for pure water) that, in turn, further diminishes the probability of solid contact and wear at the interface.
The fluid forces exhibit 'jumps' at the switching instances, which can be explained by the instantaneous change in the boundary conditions of the EHL model: during the upstroke, a large pressure differential p 23 is maintained at the piston-cylinder interface that goes to zero when the piston flaps open in the downstroke. As with the case of the fluid column mass, these artificial discontinuities will be removed in future iterations of the model by accounting for the dynamical behavior of the piston flaps. While most time histories presented in this paper, including that of the shear force in Fig. 6(c), are roughly periodic with the sinusoidal steady state having been reached by the second cycle, the time history of the normal (bearing) force in Fig. 7(c) is non-periodic: this partly arises from the fact that the pendular mass changes significantly between the up-and downstrokes, as captured in Equations (35) and (39), which in turn affects the oscillations in the x-direction. Since the shifting pendular frequency is not in phase with the heave motion, the result is non-periodic.
Another observation arising from Fig. 7(c) is the relatively small amplitude of the fluid bearing force F n , which is defined in Equation (28) as the component of the radial force acting along the x-direction. Since the fluid film is assumed to exist within the entire piston-cylinder clearance, the magnitude of the normal force,   which depends on the small difference in pressure between the positive and negative x-directions, results in a correspondingly small normal force. For example, when x > 0, integrating the pressure distribution over the 'positive' piston face (Àp/2 < q < p/2) and taking the force component along the x-direction yields a net normal force that pushes the piston in the negative x-direction; doing the same for the pressure distribution along the 'negative' piston face yields a comparable force pushing the piston in the positive x-direction. The difference between these forces gives rise to a net normal force that has a comparatively smaller magnitude and opposes the direction of motion along the x-direction.
The piston displacement and velocity are plotted again in Fig. 8(a) to show that the higher-frequency vibrations in the velocity give rise to corresponding oscillations in the acceleration time history plotted in Fig. 8(b). Their amplitude is maximum at the switching instances as would be expected. The pumping force derived in Equation (18) comprises pressure-difference and inertial contributions and is nonzero only during the upstroke when the piston flaps remain closed. Introducing a nonzero wave damping coefficient B dampens some of the higher-frequency vibration content and allows for the acceleration to become almost zero at the top-and bottom-dead-centers of the piston oscillation; the same can be achieved in a more realistic manner by introducing the piston flap dynamics into the model. Consequently, the time histories of piston velocity and acceleration, along with the relevant forces such as the pumping force, would become smoother around switching instances. In the absence of the above, we utilize smoothing of the fluid column mass calculation to remove the detrimental effect of these non-physical higher-frequency vibrations from the simulation results (Fig. 9, below).
Furthermore, it can be observed that the magnitude of the pumping force is relatively small compared to the net force on the buoy (~30 kN compared to~1 MN, respectively); however, while the order of magnitude of the net force on the buoy is expected to remain more or less constant, the pumping force can change significantly depending on the pumped volume. Specifically, it is straightforward to derive how scaling the piston radius by a factor X will scale the pumped volume by a factor X 2 and how, for example, a piston of twice the nominal radius will pump four times the working fluid volume increasing the pumping force to~100 kN. For a multi-piston pump designed to employ three pistons with radii of R 1 ¼ 5 cm, R 2 ¼ 2R 1 ¼ 10 cm and R 3 ¼ 2R 2 ¼ 20 cm simultaneously when very large waves come into the WEC, the total working fluid column volume would be~PX 2 ¼ 0.5 2 þ 1 2 þ 2 2 ¼ 5.25 times the nominal volume and, correspondingly, the total pumping force would be~525 kN with potential peaks at twice that magnitude during switching instances. Such forces would be sufficient to fully submerge the buoy, would alter the system's resonance and should, therefore, be accounted for in the MPP design; the latter lies outside the scope of the present work.
An exponential growth/decay function was introduced into the formulation of the fluid column mass as shown in Equations (19) and (20). Its effect is shown in Fig. 9: the switching between the upstroke and downstroke is done passively via the criterion of the piston velocity sign as shown in Fig. 6(aeb) and reproduced in Fig. 9(a) below.
The higher-frequency vibrations result in the rapid switching observed around seconds 15, 25, 25 and 45, and induce nonphysical vibrations as discussed previously. By using the exponential growth/decay function, switching is assumed to take some amount of time dt that depends on the value of the growth/decay   rate and the nominal fluid column mass, rather than occurring instantaneously between iterations iÀ1 and i. For example, using a growth/decay rate of G¼50 kg/s and nominal fluid column mass m f,nom z 3,500 kg, the fluid column mass goes from 0 to m f,nom or vice versa within This is shown in Fig. 9(b) where the behavior of m f has been smoothed relative to the switching. In future work, we plan to investigate the dynamics of the piston flaps and remove artificial smoothing altogether.
The net instantaneous flow-rate of a single-piston pump is plotted in Fig. 10(a). As per the model formulation, the working fluid gets pumped to the upper reservoir only during the upstroke and the flow-rate becomes zero otherwise. The net flow-rate accounts for the leakage, given in Equation (25) and plotted in Fig. 10(b), which is~1% of the gross value for the pistoncylinder clearance of 400 mm used in the simulations. This is given by Fig. 10(c) is a plot of the leakage velocity given in Equation (26) and provides an indication of the loss of pressure from the pistoncylinder separation. Leakage through the piston-cylinder clearance occurs during the upstroke when the positive leakage flow-rate is subtracted from A c v z to yield Q net . Volumetric leakage is very sensitive to viscosity. When using pure water, for example, its lower viscosity results in significant leakage for the same piston-cylinder separation of 400 mm; in order to reduce this to the levels achieved with an HFC lubricant, the separation must be reduced to 100 mm, which significantly increases the probability of solid contact (as shown in Figs. 6(d) and 7(d)).
The transfer of working fluid from the lower to the upper reservoir due to the pumping action of the reciprocating piston can be quantified as the change in pressures p 1 and p 4 as shown in Fig. 11(aeb). The instantaneous hydraulic head in each container can be calculated as the ratio of these pressures to the product of the working fluid density and the gravitational constant, i.e. L U ¼ p 1 /rg and L L ¼ p 4 /rg, respectively. As would be expected, the pressure (and head) increases in the upper reservoir and decreases in the lower by an amount corresponding to the volume of pumped working fluid. The maximum hydraulic head available for the turbine is plotted in Fig. 11(c) and includes the height difference between points 1 and 4 (refer to Fig. 2); this increases by roughly 5 mm per cycle.
The instantaneous pressures p 2 and p 3 above and below the piston flaps are plotted in Fig. 12. A large pressure difference p 23 is developed almost instantaneously at the closing of the piston flaps, while the two pressures converge to a common value during the downstroke when the piston flaps open. The relevant formulations for these pressures are given in equation set (24) and Equation (27), while the effect of this instantaneous switching was discussed previously within the context of the boundary conditions used in the EHL model. Correspondingly, nonphysical oscillations close to the switching instances can be removed in future models by accounting for the piston flap dynamics.

EHL pressure distribution at the piston-cylinder interface
Referring to the piston displacement time histories in the vertical and lateral directions, given in Figs. 5(b) and 7(a), respectively, we examine two instances of EHL pressure distributions corresponding to local maxima and minima of the lateral displacement.
The fluid film pressure distribution at the piston-cylinder interface exhibits localized pressure buildup, albeit of very small Fig. 11. Pressures in the upper (a) and lower reservoirs (b) and the (maximum) hydraulic head (c).  magnitude (<150 Pa), as shown, for example, in the downstroke instance plotted in Fig. 13(a) and seen in similar numerical experiments [22]. Here, the pressure differential p 23 is zero and p 2 and p 3 have converged to a value around 1072.35 kPa. By definition, the minimum film thickness and maximum film pressure will occur along the x-axis as this was defined in Fig. 3: for positive x, at q ¼ 0 as is the case for the plotted instance, and for negative x at q ¼ ±p (not shown).
During the upstroke, the inlet pressure is equal to p 2 and the outlet pressure, which never becomes smaller than the vapor pressure of water thereby avoiding cavitation, is equal to p 3 as per the boundary conditions of equation set (23). A pressure spike cannot be resolved in the plot of Fig. 13(b) due to the large pressure difference between the inlet and outlet conditions; this could presumably be achieved by using a much finer mesh for z and q at immense computational cost.
As stated previously, the pressure buildup in the working fluid appears to be sufficient to provide bearing support for the normal forces developed during the lateral piston motion without leading to solid contact: this is the case when using HFC, but not when using pure water, and suggests that HFCs could successfully reduce wear. This behavior could be much more complex if viscous, thermal, transient and other effects were accounted for at the piston-cylinder interface. These will be investigated in subsequent models.

Power generation and loss
Power generation and loss in system components can be calculated over the duration of simulation time using the general formula P ¼ F$v. Hence, the net power of the buoy due to the buoyancy (3), drag (6) and excitation forces (7) can be found as follows: while the power corresponding to the spring (stiffness) and damping terms of the connecting rod is with the work or change in energy e defined as the integral of power over time e being conserved for the stiffness term but lost in the damping term. Additional power losses will occur at the pistoncylinder interface due to EHL (or solid contact and friction) forces and can be calculated as The net power of the buoy fluctuates during the wave period with the power generated and lost, respectively, due to the buoyancy, drag and excitation components of the net buoy force as shown in Fig. 14(a). The current formulation of the excitation power affects the periodicity of the net power, whereas the buoyancy and drag powers show the same period as the excitation. In reality, significant buoy power will be lost due to the wave damping term of the excitation force but this is neglected in the present model The wave energy flux for the waves under consideration (with significant wave height H m,0 ¼ H w ¼ 4 m) can be calculated per meter of wave-front as P wave y0:5H m;0 2 T w ¼ 80 kW=m. For the width of the buoys W b ¼ ffiffiffiffiffi ffi A b p ¼ 7 m and the wave period T w ¼ 10 s, the total energy in the wave can be calculated as E wave ¼ P wave W b T w ¼ 5600 kJ. Out of this energy, only~365 kJ are transferred to a single buoy (assuming zero wave damping), calculated as the integral of the net buoy power over one cycle: hence, the hydrodynamic efficiency of a single buoy is small at 6.5%. The second buoy will extract energy from a smaller wave, the third one from an even smaller one, and so on as the wave energy is gradually absorbed by the device that will employ an array of adaptable pumps to maximize energy extraction. Returning to the considerations for a single buoy, in the limiting case where the total pumping force (including inertial and pressure difference contributions) will be equal to the total force that a buoy can support (~1 MN), then the total energy that could be extracted by each SPP would be roughly equal to the pumping efficiency times the buoy energy, i.e. 98.5% Â 365 kJ ¼ 360 kJ. Since it is undesirable to have a situation where the pumping force is close to the net force on the buoy, a larger number of pumps adapted to the incoming wave profile will be necessary to ensure that the buoy will always remain semi-submerged and operating optimally at resonance. Doing so will also require the use of robust control algorithms that will be investigated in future work together with physics-based characterizations of interactions between buoys.
Energy is transferred from the buoy to the piston via the connecting rod and from the lower to the upper reservoir via the working fluid. Even in the presence of buoy-piston damping, the energy in the rod is mostly conserved as shown in the power time history of Fig. 14(bec) since P C,rod is small compared to P K,rod . The power lost at the piston-cylinder interface due to EHL, plotted in Fig. 14(c), is similarly small. This power, which is zero around the switching instances, and the corresponding energy loss would be dissipated primarily as heat. Hence, an improved EHL formulation accounting for thermal effects becomes important for future work.
Pumping power is calculated as with the pumping force having been defined in Equation (18). The pumping energy (integral of the pumping power over time) is transferred via the working fluid and stored as potential energy. The potential power and energy are given in the following set of equations and are functions of the (maximum) dynamically varying hydraulic head H, shown schematically in Fig. 1(a), and the volume of working fluid pumped (Q net ¼ _ V net ), with the net flow-rate having been defined in Equation (48).
The pumping and potential powers are close in amplitude, as shown in Fig. 15, with the deviation being attributable to the pistoncylinder interface (EHL) power loss plotted in Fig. 14(c) and the losses due to the transfer of momentum of the moving fluid column to the stationary fluid in the lower reservoir. The deviation between the corresponding pumping and potential energies (areas under the curves in Fig. 15) yields the efficiency of the pumping system which is 98.5% for the current configuration with piston-cylinder separation s ¼ 400 mm, but without accounting for losses due to wall friction. To do so, we can employ the Darcy-Weisbach friction factor quantifying hydraulic loss due to wall friction that can be solved for in the laminar flow regime as f ,885 for the largest (four times larger diameter), assuming a velocity of 1 m/s, as shown in Fig. 6(b). Correspondingly, the Darcy-Weisbach factor will vary between 0.044 for laminar flow and 0.036 for turbulent flow under the assumption of a (rather high) roughness amplitude ε ¼ 50 mm. The ratio of head loss to the length of the hydraulic cylinder will be and its value is 0.022 for the laminar case and 0.0046 for the turbulent case for a moderate velocity of v z ¼ 1 m/s. In general, the influence of wall friction can be to reduce the mechanical efficiency due to hydraulic losses by as little as~0.5% or as much as~2.2% for the examined velocity.
It should be noted that these calculations assume that the cylinders span the entire height of the working fluid column (~200 m); however, this need not necessarily be the case, since the cylinders could be limited to the length of the expected maximum stroke (~12 m), and the remaining height could utilize ducts of larger diameters to transfer the fluid from the lower reservoir and into the upper one. This would result in a lower head loss Dh for the expanded duct sections due to their increased diameter (and Reynolds number) and reduced length, even though minor hydraulic losses would also be introduced due to flow expansions and constrictions. The general issue of hydraulic losses will be addressed in future pumping system designs and corresponding models.
We performed two additional simulations to investigate the effect of wave frequency on pumping behavior: one for short period waves (T ¼ 4 s, l ¼ 25 m and H ¼ 1 m) and one for long period waves (T ¼ 20 s, l ¼ 625 m and H ¼ 6 m). It should be noted that, while the connection between period and wave length is given in linear wave theory (l ¼ gT 2 /2p), there is no clear relationship between the period and the wave height; hence, we arbitrarily chose the heights of 1 m and 6 m for these two simulations, respectively. We found that the energy extracted per cycle, as a representation of the hydraulic head added to the upper reservoir, scales roughly with the wave height: a 1-to-6 ratio of 25 kJ per cycle for 1 m-tall waves versus 157 kJ per cycle for 6 cm-tall waves. The mechanical efficiency remains comparable to that of the original case we studied at 98.6% and 97.4% for the short and long wave periods, respectively. Piston velocity was within the original envelope of ±1 m/s for the two simulations and, as a result, the peak values of the leakage flow rate and leakage velocity remained constant at around 0.5 L/s and 2 m/s, respectively. Hence, the volumetric efficiency also remains constant for longer wave periods, but will increase with increasing net flow rate for taller waves. These findings suggest that the system will work over the expected range of wave frequencies.
The efficiency calculations were validated by experimental measurements on an experimental prototype of the SPP where pure water at 20 C was used as the lubricant [37,38]; these are summarized in Appendix 2. It was observed that the volumetric leakage occurring at the piston-cylinder separation is the main determinant of the pumping efficiency eitself a function of the mechanical and volumetric efficiencies of the pumpe, which decreases to below 80% for separations larger than 200 mm. The same results were validated with finite element simulations (COMSOL Multiphysics software) [23] and point to target piston-cylinder separations around 100 mm when pure water is used as a lubricant.
As stated previously, using a more viscous lubricant such as HFC allows for increasing separations while maintaining high efficiencies. It is expected that friction and lubrication at the valve interfaces will further reduce this efficiency, something that will be explored in future work.
The maximum instantaneous power that can be generated per cycle corresponds to the maximum hydraulic head available for the turbine and is found at the peak of the second upstroke (selected so as to avoid initial transients) at 30 kW; correspondingly, the potential energy generated per cycle is 99 kJ for a single piston pump. Multiple pistons will scale up energy production for a single pump unit and will be able to adapt to the energy content of the wave, while this behavior is expected to improve accordingly for the MP 2 PTO WEC. Overall, the use of multiple pump units is designed to incrementally extract almost all of the energy available in an incident wave and can transfer this energy to the stored working fluid with a high efficiency.

Conclusions
The mechanical design of the MP 2 PTO WEC, together with a tribological model of a single-piston pump, have been presented for a novel renewable energy harvester termed the Ocean Grazer. The model accounts for the pump's dynamics, hydraulics and lubricated contact under a number of assumptions and gives estimates of the maximum potential energy gained per stroke. Within the full-scale Ocean Grazer system utilizing many such pumping units, this energy can be stored without losses and used, for example, to generate electricity on demand. EHL has been shown to dominate at the piston-cylinder interface, suggesting that water-based hydraulic fluids can potentially be used as a lubricants for this application with mechanical efficiencies close to 99%, and cavitation was avoided in the proposed pump design.
Future work, as summarized below, is necessary to relax the assumptions used in the present model: Investigation of wave damping and drag as a function of buoy design; EHL for pressure-dependent lubricant viscosity, thermal and transient effects; Improved formulation of piston flap and check valve dynamics to remove discontinuities; Calculation of the drag force acting on the piston while sinking into the water column; Inclusion of additional degrees of freedom for the buoy (e.g. downstream translation, roll) and piston (pitch), and accounting for deformations and the accurate characterization of damping in the connecting rod/cable; Investigation of tribology at secondary interfaces such as the rod/cable seals; Formulation of realistic wave displacement profiles accounting for the loss of energy content due to energy extraction for use with multiple buoys in a grid (floater blanket).
The present model will be used as a building block in the development of a multi-piston pump model coupled to variableload control. A scale prototype of the multi-piston pump has been recently constructed and was used to validate model predictions: indeed, the reported mechanical efficiency of the system is close to 99%, with the pumping efficiency (h pump ¼ h mech Âh vol ) being dominated by volumetric losses at critical interfaces. In turn, the resulting model can be used to optimize component and system design to increase power generation and reduce losses. In addition to parallel research activities in the hydrodynamic design of the buoy grid (floater blanket), and the floating structure, this work is fundamental in the realization of the Ocean Grazer and is a necessary tool in the performance of feasibility studies focusing on energetics, cost analyses, coupling to existing power grids, etc. : The difference in the radial distance from the origin for the cylinder and piston respectively yields an analytical expression for the film thickness where w r denotes the deformation of the solid surfaces and sin4z4 and cos4z1 for a very small tilting angle, while L COM ÀH p z L COM and Àp q p. The solid surfaces of the piston and cylinder will deform even at the lightest loads allowing for the fluid film to provide the essential bearing load capacity [31,33]. Assuming only elastic deformation, w r must satisfy the following equation: where the z-coordinate lies in the direction of z and the combined elastic modulus is defined as : (A. 5) Equation (A. 4) assumes smooth elastic surfaces. Roughness could potentially be incorporated in the prediction of the deformation utilizing multi-asperity [34,39e44] or fractal representations of surface roughness [45,46]. Similarly, finite element analysis could be used to predict the displacement fields of the piston and cylinder surfaces [22]. Given that the piston-cylinder separation will be larger in the proposed single-piston pump than in other applications such as reciprocal engines [20] and ring-less compressors [22], deformations are expected to remain within the elastic regime. Hence, Equation (22) is deemed sufficient as a first approximation of deformation in the EHL problem.
The Gauss-Seidel iteration method is used to solve the EHL problem [35]. The discretized form of the Reynolds Equation (22) becomeswhere the partial derivatives are defined as follows: vp vz methods, as well as the experimental and simulation parameters used, can be found in the relevant conference paper [37], but representative results are provided here for completeness. It should be noted that pure water at 20 C was used as a lubricant in, both, the experiments and simulations while the lubrication model was simplified with the use of hydrodynamic approximations for the normal and shear forces instead of the full EHL solutions discussed in the present paper.    A1(a) shows the kinematics of the SPP prototype: the displacement was measured with a displacement sensor (CELESCO SP1-50) and the velocity was calculated by numerically differentiating the displacement time history. The pumping force in Fig. A1(b) was measured with a force transducer (TEDEA HUN-TLEIGH 615) installed at the cable connecting the piston to the actuating mechanical arm that was used to simulate sinusoidal wave motion, while the power production in Fig. A1(c) was calculated as the product of the pumping force and the velocity. The same results were predicted with the simplified SPP dynamical model whose parameters were modified to match the experimental conditions, as can be seen in Fig. A2. From the comparison of the two, it becomes clear that the oscillations observed at the switching instances of the simulated results are an artefact that must be corrected for in the dynamical contact models, as discussed in the body of the present paper: such oscillations might occur at the beginning of the downstroke, as can be seen in the experimental results, but their amplitude is small compared to that predicted in the simulations.