Fully Nonlinear Time Domain Analysis for Hydrodynamic Performance of an Oscillating Wave Surge Converter

The hydrodynamic behaviour of an oscillating wave surge converter (OWSC) in large motion excited by nonlinear waves is investigated. The mechanism through which the wave energy is absorbed in the nonlinear system is analysed. The mathematical model used is based on the velocity potential theory together with the fully nonlinear boundary conditions on the moving body surface and deforming free surface. The problem is solved by the boundary element method. Numerical results are obtained to show how to adjust the mechanical properties of the OWSC to achieve the best efficiency in a given wave, together with the nonlinear effect of the wave height. Numerical results are also provided to show the behaviour of a given OWSC in waves of different frequencies and different heights.


Introduction
With the extensive use of fossil energy, it becomes evident that this is not sustainable as the energy resource can be exhausted. Furthermore the excessive use of this kind of energy has caused environment pollution and climate change. Therefore there has been an increasing interest in developing technology to efficiently make use of renewable energy, in the form of water wave, hydro, wind and solar powers, etc. Indeed some of technologies are already widely in use, and had supplied an estimated 22.1% of global electricity consumption by the end of 2013, among which hydropower provided about 16.4%, and wind energy produced around 2.9%. In contrast, use of wave energy is less developed while it has a great potential.
Most wave energy exits in deep seas, and thus the earlier wave energy converters were usually designed for deeper waters (Salter, 1974. While during recent years, the wave energy in nearshore has received increasing attentions, which is a result of the new understanding of the nature of wave energy. It is common that the nearshore wave power resource is dramatically less than that offshore. However this statement is usually based on the omni-directional or gross wave power resource, which is particularly advantageous for the devices which are not sensitive to the direction of wave propagation such as isolated axi-symmetric heaving buoys. However for the devices which respond differently to waves from different directions, "the exploitable wave energy resource", which is defined as the mean value of the directionally resolved incident wave power, where the threshold power level is not allowed to be more than four times the mean value, is a more realistic measure . Based on such a new definition, there is only 10-20% energy losses from the offshore to the nearshore in many sea sites . From engineering and economy perspectives, the devices located in the nearshore have a number of advantages. Through the depth induced wave breaking, the large although infrequent extreme waves, which may be a threat to the survival of wave energy devices, can be filtered (Henry et al., 2014). In addition, cost of the cables which are used to bring power to the shore and power losses in the cables will be lower due to closer distance to the coast. The installation and maintenance costs in nearshore are likely to be reduced due to the shorter weather windows, the closer location to the shoreline and shallow water depth (Folley et al., 2007a). Thus the wave energy converters designed for shallower water gradually gain momentum during the last decade. This type of devices are often classified as oscillating wave surge converters (OWSCs), as energy is mainly extracted from the horizontal or surging motion of water particles in ocean waves (Whittaker and Folley, 2012). One of the most typical OWSCs is Oyster. This device is a buoyant and hinged flap which is attached to the seabed and has the top edge piercing through the water surface. It is usually used in area of water depth of between 10 and 15 m, which is often referred to as nearshore. Two full-scale prototypes of Oyster were installed in summers of 2009 and 2011 at the European Marine Energy Centre's test site (EMEC), Orkney, Scotland. Apart from Oyster, WaveRoller (Folley et al., 2007b) with a flap completely submerged is another typical representative of OWSCs. The previous two devices typically comprise a single flap spanning their full width. Recently, a new concept 'Modular Flap', formed by splitting the flap into multiple modules, is receiving increasing attention. Sarkar et al. (2016) made comparison of the performances between a rigid flap and a modular system with six identical modules. The combined length of the latter was the same as that of the former. It was found that the module system had the potential to gain more energy due to its multiple natural frequencies at which resonance could occur. Wilkinson et al. (2017) showed through both the numerical and experimental methods that the modular device could absorb power more smoothly.
Over the past decade, OWSCs have been investigated by many authors though various analytical and numerical methods. Folley et al. (2005Folley et al. ( , 2007a investigated the effect of water depth on the power capture efficiency and surge wave force of an OWSC through the combined method of analytical and numerical results. The capture efficiency and surge wave force were obtained through the equations with hydrodynamic coefficients obtained from the linear potential flow theory. Similar work includes those by Folley et al. (2007b), Van't Hoff (2009), andFolley (2012). In addition, a three dimensional semi-analytical model based on the potential flow theory was developed by Renzi and Dias (2012a, 2012b to investigate the behaviour of an OWSC. Recent years, generally CFD techniques have been increasingly used. Among them, Smoother particle hydrodynamics (SPH) method was used by Rafiee et al. (2013) to investigate the wave loading on an OWSC, and by Yeylaghi et al. (2016) to investigate the interaction between wave and OWSC. Another important CFD technique is to combine the finite volume method (FVM) in solving RANS function, with some surface tracking methods, such as volume of fluid (VOF) or level-set (LS) technique for the free surface. In this context, viscous effects and slamming effects on an OWSC were investigated by Wei et al. (2015Wei et al. ( , 2016 though the Fluent Software. Similar method was also adopted by Schmitt and Elsaesser (2015) in simulating the motion of an OWSC in significant waves through an open source toolbox OpenFOAM.
Within the linear theory, it is generally accepted that the wave energy absorbed by WECs from waves very much depends on the difference between the wave and body frequencies and the damping from the power take-off system. Specifically, the wave energy converter must first be resonated at the design wave frequency. Simultaneously, the radiation damping of the WEC should be matched with the externally imposed damping. Mei (2012) named it as "impedance matching", and found that the power output could be maximized when "impedance matching" was realized. Therefore achieving resonances is recognized as one of the best options to improve the performance of a WEC. For instance, a coupled resonance was considered by Evans and Porter (2012) through inserting a water tank having separate mass /spring /damper characteristics into a WEC. The multiple resonances of system were exploited by Crowley et al. (2013) through inserting a mechanical system of compound pendulums inside the hollow cylinder to provide a broad-banded response. Other similar work includes those by Eriksson et al. (2005) and Renzi and Dias (2012a, 2012b. Within the potential theory, the linear frequency domain method, which is most commonly used above, is very effective for the efficiency evaluation of WECs in small regular waves. However, in large waves, the neglected nonlinear effects may lead to the significant deviation of the efficiency and performance evaluation. In such a case the fully nonlinear time domain method may be needed. In the present paper, we shall focus on the oscillating wave surge converter undergoing large amplitude motion. The governing Laplace equation will be solved by the boundary element method (BEM). The nonlinear free surface boundary conditions will be satisfied through time stepping method. In the following sections, we shall first introduce briefly the OWSC system, with the definitions and equations for the absorbed power and efficiency. This is followed by the detailed description and the methodology for hydrodynamic analysis. Results are then provided to assess the performance of the device in the nonlinear periodic waves. Efficiency analysis is undertaken in a given incoming wave while varying the natural frequency and the mechanical damping of the system. Analysis is also made to access the nonlinear effect through the incoming wave amplitude. All these aim to provide a better understanding of the mechanism of wave energy conversion.
2 Methodology for efficiency prediction 2.1 Definition of wave energy absorbing efficiency of an oscillating wave surge converter Fig. 1. The sketch of an oscillating wave surge converter From a practical point of view, OWSC is a three dimensional device. However two dimensional analysis can still shed some insight into the physics of the problem. Fig. 1 gives a sketch of a two dimensional OWSC (Henry et al., 2014), in which a flap hinged at distance h from the calm water surface, is allowed to rotate about a horizontal axis perpendicular to the direction of wave propagation. The Cartesian coordinate system O xy  with the origin O fixed at the hinge centre is defined so that x axis is parallel to the flat seabed with depth d and y axis points vertically upwards. The angular displacement  is zero when the flap is in the vertical position and is positive in the anticlockwise direction.  is the angular velocity and therefore / d dt . In the present paper, the distance h , the acceleration due to gravity g , and the water density  are used for dimensionalisation. The flap will be in oscillatory rotation under the excitation of the wave, and its hydrodynamic behaviour is to be investigated through the time stepping method. It is expected that the rotation of the flap is resisted by the forces exerted by the generator. Thus the equation of motion of the flap can be written as where dot denotes temporal derivative, I is the rotational inertia of the flap about the hinge, pto a and pto c respectively imply the inertia and elastic characteristic of the generator, while pto b implies the mechanical damping due to power take off system or generator.
which shows that pto b is the equivalent energy extraction rate.
To consider efficiency of the device, we may consider the energy propagation in the incident wave. If the wave energy E confined within 0 x  , where  is the wavelength, we then have (Wehausen and Laitone, 1960) The first term is in fact the energy leaving the domain and the second term is the energy moving into the domain and is the one to be absorbed. We notice that for the progressing wave, in which cT   has been used and  is the elevation of the free surface. The efficiency R of the system can be defined as the ratio of the power extracted by the device to the power in the incident wave, or 2.2 Hydrodynamic analysis for fully nonlinear interaction between waves and converter For the wave problem, water can be assumed to be inviscid as the viscous effect becomes important only after many wave periods or wavelengths (Lighthill, 1978). When the fluid is further assumed to be incompressible and the flow to be irrotational, the velocity potential whose gradient is equal to the fluid velocity can be introduced and it satisfies the Laplace equation. Let the incident potential be denoted by I  and the incident wave elevation by I  . Eq. (6) becomes In addition to the incident wave, the moment M in Eqs. (1) and (2) The dynamic and kinematic boundary conditions on the free surface F S in the Lagrangian form can respectively be written as At the far field, the disturbed wave will propagate to infinity and there should be no reflection.
Numerically, the disturbed wave will be absorbed and then diminishes when x .
In such a case, the potential will tend to the incident potential, or To implement in the numerical simulation, this condition is imposed at a truncated vertical boundary at a sufficiently large distance x . On the free surface near the truncated boundary, a numerical damping zone, is introduced in the dynamic and kinematic free boundary conditions in the following mixed-Eulerian-Lagrangian (MEL) form (Contento et al., 2001, Aliabadi et al., 2013) includes the variation of the potential due to the change of the free surface elevation.
We notice that the last terms in both Eqs. (15) and (16) correspond to the difference between the total potential (total wave elevation) and the incident potential (incident wave elevation). This is to ensure only the disturbed potential, including the diffracted and radiated potentials, is absorbed, but not the incident potential I  . The damping coefficient () x  in these two equations is chosen such a way to ensure a smooth transition from Eqs. (12) and (13) to (15) and (16) where  and  mean the frequency and nonlinear wavelength of the incident wave respectively. The strength and the length of the damping zone are respectively controlled by two coefficients  and  . Based on the numerical tests in the presented problem, it is found that values of  and  can both be taken as 1.0. For initial condition at 0 t  in the mathematical modelling, the flap can be assumed to be put in the incident wave suddenly, and the potential on the free surface and the wave elevation can be taken as The mathematical problem is solved using the boundary element method. Through Green's identity, the Laplace equation in the fluid domain can be converted into an integral equation over its closed boundary S . Within each element the potential and its normal derivative are assumed to vary linearly. Boundary conditions are used, and the unknowns are then found from the solution of the matrix equation (Lu et al., 2000, Sun and. Once the solution of the above problem is found, the fluid moment may be obtained from the integration of the local moment created by the pressure over the flap surface. We have where m is the mass of the flap and c y is the vertical coordinate of the gravity centre of the flap when 0   . We note that in the equation, even when the velocity potential  has been obtained at each time step numerically, t  is still not explicitly known. An effective method to resolve this is to treat t  as another unknown function which satisfies the Laplace equation in fluid domain. On the body surface, the normal derivative of t  can be written as (Wu, 1998) We note that in the equation the acceleration   is not yet known before Eq. (1) is solved, which in turn depends on M in Eq. (21) and   in the Eq. (22). To decouple this nonlinear mutual dependence of the flap acceleration and fluid loading, we adopt the method of Wu andEatock Taylor (1996, 2003). In particular, we define 1  and 2  as and their numerical calculations can be performed based on the procedure in Xu and Wu (2013). Based on the Bernoulli equation, together with zero pressure on the free surface, we have At 0 t  , the free surface and the potential on the free surface are prescribed based on the incident wave, together with the initial value of  and  . The solution of the discretized form of Eqs. (19) and (31) then gives   , from which  and  can be updated though the Runge Kutta process. At the same time, Eqs. (12) and (13) are used to update the potential and free surface. The calculation then moves to the next time and it continues until the desired time step.
On the free surface the discrete element nodes may cluster or scatter in the process of updating, and thus the mesh needs to be re-generated after every few time steps. Here we choose the 5-point 4 th order Lagrangian interpolation method to redistribute the nodes (Abramowitz and Stegun, 1964) is the Lagrangian interpolation polynomial. We adopt the length measured along the free surface as variable e , with its origin at the intersection between the free surface and body surface. xy at e .
The present methodology can deal with any given incoming wave. The tenth order and fifth order Stokes waves may give similar overall results in the present cases, such as total force and motion. When the detailed results such as higher order force components are needed, the higher order Stokes wave is more appropriate. Thus we choose the 10 th order periodic Stokes wave through Fourier approximation method (Fenton, 1988) as an example. The wave elevation and potential can be written as  (Fenton, 1988). All the parameters above are nondimensionalized based on the manner described before Eq. (1).

Linear frequency domain method for fluid flow
When the wave amplitude and motion amplitude are small, the problem can be simplified as a linear one. We may follow the analysis in Section 7.9.2 of (Mei, 1983).
The last two terms are related to the displacement  and therefore together they constitute the restoring moment. The boundary conditions will be imposed at the mean position of each surface and only linear terms will be retained. Correspondingly there will be only the first term of 1 n  in Eq. (3). Also only the linear term in Eq. (35) where b m is the mass of the water displaced by the flap at the vertical position, b y is the vertical coordinate of the buoyancy centre and sin  has been used due to linearization. Eq. (1) can then be transformed as For the linear problem and a symmetric body, the Haskind relationship (Mei, 1983) For the linear wave, Eq. (6) gives 2 g 1 2 Thus the efficiency in Eq. (7) becomes To find the added inertia which is the same as that of Eq. (9. 10) of Section 7.9.2 of (Mei, 1983). Substituting Eqs. (49) and (56) into (7), we notice that the highest efficiency is 50%. This is a well known result for a symmetric device based on the linear theory and will be assessed below when the nonlinear theory is used.

Convergence study and validation
The computational domain is truncated at   Fig. 4. The tip of the jet was cut regularly (Sun et al., 2015). The results are compared with the experimental data of Henry et al. (2014) and those from their numerical simulation based on the finite volume approach through Fluent Software. Dimensional parameters are used in the figure to be consistent with their results. It can be seen from the figures that all the results are in good agreement in general, while the present results are even closer to the experimental data. This validates the present methodology and procedure.  , which is small enough that nonlinear characteristics will not have major effect on the results. Fig. 5 gives the motion amplitude 0  and the rate of power extracted P E against n  at different pto b . P E is calculated from Eq. (2) when the results become virtually periodic, for example at 20 tT  , and the number of periods m are taken as 20 or more. 0  is measured from the mean peak value during the m periods.
According to Eq. (53), n  can be adjusted by changing pto a or pto c , and the effects will be similar. In the current case, pto c is adjusted to vary n Fig. 5 give the results for  =0, 0.5, 1.0 and 1.5 respectively.
As expected at 0   , or when the mechanical damping is zero, the motion amplitude 0  is largest. However, this large motion does not lead to any energy extraction by the flap. In fact we have P 0 E  at 0   based on the Eq. (2). This means that all the energy has been returned to wave. To absorb energy,  must be larger than 0. It can be seen when 0.5   , the motion amplitude becomes smaller due to the additional damping. However, exactly due to the reduced motion, some of the energy has been transferred to the generator, which can be seen through the curve of P E in Fig. 5(b). Based on the linear theory, when pto b is given by Eq. (54), P E will reach its optimum at each n  . This can be seen from the curve of P H are provided in Fig. 6.
All the curves for motion amplitude in Fig. 6 Thus based on Eq. (6), the efficiency curves in Fig. 6(c) follow the exact shapes of the corresponding P E curves in Fig. 6 ,0.2,0.5 in Figs. 6(b,c) are different. Based on the discussion after Eq. (56), the highest possible efficiency from the linear theory is 50%. This can be observed from Fig. 6  In the above case, the frequencies of incident waves are all larger than natural frequency n  and there is no peak in Fig. 7. We now adjust the natural frequency n