Near wake hydrodynamics and structural design of a single foil cycloidal rotor in regular waves

We present a hydrodynamic and structural model to design a single foil wave cycloidal rotor in regular waves. The hydrodynamic part considers potential flow and represents the foil as a point vortex. The effect of the point vortices left on the wake of the foil and a correction for flow separation are considered. The structural part utilises beam theory to compute the bending moments and stresses on the foil of the cyclorotor. The validity of the hydrodynamic model is explored in attached and vortical flow conditions with the aid of CFD. Results show that the hydrodynamic model estimates the mean loading on the foil within 15% for attached flow conditions, whilst it underpredicts the loads in vortical flow conditions. Furthermore, large excursionsfromthemeanloadarefoundduetovortexsheddinginthelatter.Becausetheoptimal structural operation of the rotor is in attached flow conditions, we utilise the coupled model to design a rotor that operates optimally for a range of different sea conditions. We find that with careful dimensioning of the radius and span, power extraction in regular waves can be optimised, whilst the structural penalty is kept constant at the allowable stress level.


Introduction
Wave cycloidal rotors are a novel type of wave energy converter (WEC) that have gained a rejuvenated level of attention over the past years [41,42,5,20,35,13].Although the pioneering idea of extracting wave energy through the rotational motion of a submerged foil and the use of lift forces dates back to the early 90s [25,14], recent efforts of the Atargis group in America and the LiftWEC consortium in Europe, have brought this technology closer to commercialisation.
The concept of a wave cycloidal rotor consists of a foil that rotates under a wave.The span of the foil is oriented parallel to the crest of the wave, and through rotational motion, the foil interacts with the wave particle velocity to produce lift and drag forces.Provided that the lift to drag ratio is high, the tangential component of the lift force drives the rotation of the foil and energy is extracted.This type of WECs are classified as lift based wave energy converters [28,22].Because their operation is based on lift forces, the hydrodynamic and structural challenges encountered in the marine environment are unique, but can

Principle of operation
In this section we present the concept of the single foil wave cyclorotor used in this study.The rotor is shown in figure 1.It has a curved hydrofoil connected to a central shaft.The central shaft rotates due to the motion of the foil.
The shaft is held by bearings that are embedded in triangular frames.These frames act as the support structure and are fixed to the seabed.The rotor operates in close proximity to the free surface but the foil remains submerged during operation.
The hydrofoil has a uniform cross section along the span  and rotates following the wave orbital motion.The phase of the rotation is modulated so that it is different to that of the wave.This phase difference generates an inflow velocity This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d at an angle of attack  and hence a lift force in the hydrofoil.Provided that  does not exceed the stall angle of the hydrofoil (  ), then the ∕ ratio remains high and the tangential component of  sustains the rotation of the foil.The wave cyclorotor of this study is conceptually designed to operate in the Atlantic coast of France.According to Sierra et al. [43], the range of mean energy periods   and wave heights   in this region over a range of 41 years (1958-1999) lies within 8 to 10 s and 1 to 2 m, respectively.We compare these ranges to the values shown in the wave scatter plot of figure 2. The figure shows the data corresponding to a point in the North Atlantic at the coast of France, located at 47.84 • N, 4.83 • W. The figure shows   along the horizontal axis and   along the vertical axis.The data is available from the Ifremer FTP server and contains directional spectral wave data for 10 years between 2000 and 2010 [1].In the figure, the dominant   lies within 6 to 10 s, and the dominant   lies within 1 to 2 m.Therefore, the Ifremer database is in agreement with the observations from Sierra et al. [43].Because the Ifremer database has the data available as wave energy   and we use   in our hydrodynamic computations, we covert   values to   values through   =   , where  = 0.9 for a Jonswap wave spectrum [15].From figure 2, we select the wave design 6 8 10 12 14 16 Wave period, Tp (s) This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d conditions to be a point in the region where the highest counts are found.As such, the selected wave design conditions are   = 9 s and   = 1.2 m.

Hydrodynamic model
We consider the single foil wave cycloidal rotor that was introduced in the previous section and we show its side view in figure 3. The figure shows the single foil rotor at the normalised time period ∕  = 0 and at the azimuthal position  = 0 • .The wave direction is from left to right and the wave particle motion is clockwise [12].Hence, the rotation of the foil is clockwise as well.The rotor has a radius  and a submergence depth  0 measured from the mean sea level to the central shaft.The phase angle between the foil motion induced velocity −u and the wave velocity component v is .We refer to this phase angle as the operational phase.The relative velocity and the angle of attack on the foil are w and , respectively.The force components acting on the foil are lift () and drag (), which can be decomposed into the radial () and tangential ( ) force components through .
Considering a large span hydrofoil with uniform cross section, we can assume two dimensional flow.The foil is modelled as a single point vortex moving under the free surface.As such, we describe hydrodynamic model in the following paragraphs.

Point vortex model
The complex potential of a free vortex under a free surface was derived by Wehausen and Laitone [48].This representation has been used in the literature of wave cyclorotors to model the foils as single vortices [25,41,18] or to discretise the foils into multiple lump vortices [41,18].Here, we utilise the single point vortex representation This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d developed by Emarkov and Ringwood [18] for a two foil cyclorotor, but adapt the model to a single foil rotor.The model accounts for unsteady wake effects and in this paper, is corrected for flow separation.
Noteworthy, single point vortex methods have been predominantly used to predict surface elevation due to a foil near a free surface but have not been widely used to predict the loading of foils of a cyclorotor.For this reason, in this work, we will compare the single point vortex method results to numerical results obtained with the aid of CFD RANS simulations.This will enable us to assess the validity and limitations of the single point vortex model.
Concretely, the complex potential of a point vortex under a free surface [48] is given by where Γ() in the first term of equation 1 is the circulation of the point vortex that represents the foil, and Γ() in the second term is the circulation of the point vortices left in the wake of the foil,  is a point in the complex plane denoted by  =  + , () is the position of Γ() and is defined as () =  Γ +  Γ ,  ′ () is the complex conjugate of () and denotes a mirror point vortex that imposes the impermeability condition on the free surface, such that () is the gravitational constant,  is the wave number,  is time,  is a time parameter that determines the influence of the wake vortices of circulation Γ on the velocity field at point  and  ′ () is complex conjugate of (), where () is the position of each wake vortex.
Equation 1 can be simplified by solving analytically the integral over  of the second term with the Dawson's function.As such, Ermakov and Ringwood [18] propose the following expression: ] . ( By taking the derivative of equation 2, we obtain the complex velocity from the point vortex of circulation Γ() at point , such that This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d where the indices 1 and 2 correspond to the first and second term of equation 2, respectively.In equation 3,  1 and  1 represent the induced velocity components of the point vortex of circulation Γ() located at () at point .In contrast,  2 and  2 represent the induced velocity components of the vortices in the wake of the foil all with circulation Γ() at point , during the interval  =  −  to  = , where  is one cycle of rotation of the cyclorotor and is considered to be the maximum life of the wake as observed in our CFD numerical simulations.
The explicit terms of equation 3 are included in Emarkov and Ringwood [19] for a two foil rotor.Here, because we consider a single foil rotor,  and  are the same, i.e. the point at which we evaluate the influence of Γ() is the same point where Γ() is located.As such, only the terms related to the vortices in the wake of the foil ( 2 and  2 ) are considered.We note that  2 and  2 , account for the unsteady effect that the wake vortices impose in the near flow field and the bound circulation of the foil.Although Γ is time dependent, here and similarly to previous studies [41,42,18], the circulation of the vortices in the wake is considered constant initially, and computed through the Kutta-Joukowksi theorem (Γ = ∕).Then, their time dependant influence on the bound circulation of the foil is considered through the Dawson function.Because the function is a decay function, it accounts for a reduction of the effect of each wake vortex in the velocity field at point , i.e. during the lifetime of the wake  , the further away the wake vortex is from the foil, which occurs at  =  −  , the less its influence at point  is.
Once q from equation 3 is determined, we consider also the influence of the velocity components due to the wave and the rotation of the rotor.Figure 3 shows that in a simplified form, the relative velocity on the foil w is given by the velocity triangle formed between the wave velocity component v and the velocity due to the rotation of the rotor − .
By considering also q, then the relative velocity on the foil can be defined as To compute , let us consider the rotor of radius  from figure 3. The position of the point vortex that represent the foil is given by and and where  is the rotational frequency of the rotor.Note that the induced velocity due to the rotation of the point vortex is − , therefore we consider −  and −  .
The wave velocity  and its horizontal and vertical components (  ,   ) are determined assuming deep water linear wave equations [2], such that and where  is the wave height,  is the wave period,  is the wave number,  and  denote the position of the hydrofoil, as defined by equations 5 and 6, respectively.In the numerical computations, the wave number  is computed with the dispersion relationship [31], whilst  =   and  =   .
The angle of attack () is defined as the angle between the relative velocity  and the rotational velocity of the rotor .As such,  is given by We define positive  anticlockwise as depicted in figure 3. The lift and drag force on the foil are respectively, where   is the lift coefficient,   is the drag coefficient,  is the fluid density,  is the chord lenght of the foil,  is the span and || is magnitude of the relative velocity on the foil.The tangential force on the hydrofoil is defined as and the radial force is Because  is dependent on the angular position of the point vortex (), the average tangential force  is expressed as: The total mean torque () is and finally, the mean power output is

Angle of attack oscillations (𝛼)
In figure 4, we present the angle of attack oscillations () throughout one period of revolution at Δ = 90 • for the single foil of the cyclorotor, as predicted by the single point vortex model.It can be seen, that the maximum angle of attack oscillation occurs when the foil is closest to the water surface.Furthermore, the oscillations are asymmetric, i.e. they deviate more from the mean  in the first half of the cycle and less in the second half.This is reminiscent of the oscillations that occur in cross flow turbines [4,30].Although in the example of figure 4 ,the amplitude of  is not This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d severe, it is likely that under different wave conditions the oscillations of  increase and dynamic stall -the ability of a rotor to maintain lift through severe angles of attack and in the presence of vortex shedding [39] -occurs.Hence, because dynamic stall is not accounted for in the single point vortex model, a case where the flow remains attached and a case where dynamic stall is likely to occur will be investigated in this paper with CFD simulations.This with the aim to assess the validity of the point vortex model in different scenarios.This type of assessment has not been carried out in the literature of wave cycloidal rotors.As a second novel contribution of this work, we enhance the single point vortex model by incorporating a correction for flow separation.We describe the correction in the following paragraphs.

Flow separation
To account for flow separation effects, we utilise the modified version of the Leishman-Beddoes method, as implemented by [39,10].In summary, the correction establishes a relationship between the trailing edge separation point and the normal force coefficient   on the foil.In the case of the foil of the cyclorotor,   is the radial force coefficient   , as per the notation used in equation 15.The correction is implemented through Kirchoff's theory.The coordinate of the separation point measured from the leading edge of the foil   is normalised by the chord length , such that  =   ∕.A fully attached boundary layer yields  = 1, whilst a fully detached one yields  = 0.The radial force coefficient is then computed as where  0 is the lift curve slope of the foil and  0 is the angle of zero lift.Although  0 and  0 data are not available for curved foils, it is expected that a curved foil generates zero drag as it moves along the circumference of rotation of the cyclorotor.As such, the curved foil is expected to behave as a symmetric foil in straight flight.Therefore, we consider

CFD model
The point vortex method applied in this work represents a numerically efficient approach to design wave cycloidal rotors.However, the consideration of a number of effects in cyclorotor hydrodynamics such as surface drag, wake vorticity and flow separation is of empirical nature.In order to obtain an understanding of the reliability of the point vortex method in design of the cyclorotor, we conducted CFD numerical simulations.
The simulation of wave energy converters in CFD based numerical wave tanks allows the advantage that aperiodic, non linear and viscous effects are inherently considered.With increasing availability of computational resources, these methods have increased in popularity for the investigation of specific hydrodynamic effects of wave energy converters [49].As such, the numerical model employed in the presented work corresponds to the setup described and validated in [34].
The fundamental equations which form the basis of the CFD model are the Navier Stokes and the continuity equation, such that respectively.In equation 20, the two terms on the left hand side represent the local and the convective acceleration, respectively, and ∇ is the del operator and  is the velocity vector.On the right hand side of equation 20, −Δ is the pressure gradient,  is the viscous term and   represents the body force.The formulations shown in equations (20) and ( 21) assume the fluid as incompressible, and they describe the conservation of momentum and mass in the fluid, respectively.
The solution of the equations is achieved by employing a finite volume based approach.The numerical domain is discretised by a finite number of cells and the equations are solved for pressure and velocity in each cell.As described in [34], the Reynolds averaged form of the Navier Stokes equation (RANS) can be employed for modelling of the cyclorotor, as it allows to significantly reduce the computational time.Turbulence modelling is performed through a modified version of the  -model originally presented in [33].The model is extended by the turbulence limiters This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d proposed in [29] and [16] to prevent an overproduction of turbulence in the vicinity of the free surface and near the stagnation points of the foil.
The multi phase flow problem is treated using a Volume of Fluid (VoF) approach [26].The two phases (air and water) are treated as immiscible, while local flow properties such as density or viscosity are approximated as cell averaged values using the volume fraction .This variable indicates to which degree a cells is filled with water (1 meaning only water in cell) and is treated numerically with a transport equation in analogy with (21).
In all simulations presented in this study, the free surface interface is resolved with 120 cells per wavelength and 15 cells per wave height.Stepwise coarsening of the mesh structure away from the free surface is applied as proposed in [38].A section of the mesh in the vicinity of rotor and free surface is shown in figure 5. employed in all simulations.The radius of the overset domain is 1.4.At the interface of the domains, the solution is transferred between grids by means of interpolation (cf.e.g.[32]).
In all simulations, the background domain is defined with a total length of 8, with  corresponding to the wave length of the respective simulation case approximated based on Stokes fifth order wave theory [21].The rotor axis is located at the longitudinal centre of the domain.In the lateral direction, the domain is discretised using a single cell layer, with symmetry conditions applied on each lateral wall, thus effectively providing a two dimensional setup.At the inlet, the volume fraction and velocity of the target wave are prescribed as a Dirichlet boundary condition, using the time varying values approximated again using fifth order wave theory.A sketch of the domain is shown in figure 6.
The bottom of the domain is modelled as a wall.The top boundary is defined with atmospheric pressure levels.
At the outlet, the hydrostatic pressure profile and volume fraction of a calm free surface are defined as a boundary This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d condition.A damping zone extending over a length of   = 2.5 is defined upwave of the outlet boundary to prevent reflections.The damping method presented in [37] is employed for this purpose.When employing the described setup, reflection levels of approximately 1% were found in [34], following the approach described in [46].As such, the current setup yields a minimal and satisfactory level of reflections for every tested case.

Structural model
We select a moderate strength steel for offshore applications [11], as the construction material for the foil.The mechanical properties of this type of steel are listed in table 5.In the The hydrofoil cross section is depicted in figure 7a.The foil is a curved NACA 0015 whose camberline follows the curvature of the circumference of the cyclorotor.The chord length  is measured from leading edge to trailing edge as shown in the figure.An inner skin thickness of 10 mm is considered for the walls of the foil.This is a typical thickness that is commercially available for offshore steel plates and for hollow square sections.The foil has three inner cavities, which are identified in the figure.The neutral axis is defined along the camber line of the foil and is highlighted with a black dotted line.The second moment of area  x x of any given cross section can be computed as where   denotes the distance from the neutral axis of the section to the centre of an infinitesimal area element  within the area delimited by the boundaries of the cross section.
In our example, we approximate the solution of  x x by straightening the neutral axis and approximating the shape of the straight foil with three hollow sections: triangular, square and elliptical, as shown in figure 7b.We denote the sections as section 1, 2 and 3, respectively.Each hollow section is broken down into two surfaces.A surface delimited by the outer boundary of the hollow section, and an inner surface delimited by the inner boundary of the hollow section.
In figure 7b, the base and maximum height of each of these surfaces is denoted with "b" and "h", respectively.In the figure, each base and height has two subindices.The first one denotes the outer (o) or inner boundary (d).The second one denotes the index of the hollow section (1,2,3).
Then,  x x for each hollow section is computed by subtracting the  x x of the outer and inner surfaces.Hence, the moments of inertia of sections 1 to 3 are 2 Abel Arredondo-Galeana et al.: Preprint submitted to Elsevier Page 13 of 32 This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d respectively.In equations 23, 24 and 25  o is the width of the outer sections and ℎ o is the maximum height of the outer cross sections.For all of the sections, ℎ o = 0.15.The total  x x is obtained by adding  x x,1 ,  x x,2 and  x x,3 .We note that because the  x x of all cross sections has a cubic dependence on ℎ, changes in  do not alter significantly the second moment of area.
Once the total  x x has been computed for the foil, we now focus on computing the bending moments and stresses on the foil.We model the foil as a beam and consider two benchmark types of loading: 1) uniform loading and 2) elliptical loading.The first type of loading can be promoted through large spanned hydrofoils with winglets.A large span can promote two dimensional flow, whilst winglets can prevent tip losses [24].The second type of loading is typical of elliptical planforms [27,24], as the one used in the Spitfire aircraft.This is a planform that provides a uniform lift coefficient and uniform induced angle of attack throughout its span.
The two types of loading are illustrated in figures 8. In the figure, the vertical axis shows the normalised distributed load, where    and   are the maximum amplitude of the uniform and elliptical loading, respectively.The horizontal axis shows the spanwise coordinate , where  = 0 is the origin and  =  is the span of the foil.
In order to obtain the bending moments and stresses that act on the beam, we obtain expressions for the distributed loading and forces that act on the beam as a function of .We refer to these terms as   and   , respectively.For uniform loading whilst for elliptical loading This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d The force   is the area under   .Hence, for uniform loading   =     (28) and for ellipitical loading Noteworthy, for elliptical loading the area under   has the shape of a quarter of an ellipse with a height defined by equation 27 and a base , whilst for uniform loading the shape is a rectangle of height    and base .
Equations 28 and 29 are used in the static equilibrium equations of a fixed beam to solve for the shear forces  and the bending moments .To solve for , the centroid or point of action of   is determined.For a rectangular shape (uniform loading), the centroid is at the symmetry line of the rectangle.For a quarter of an ellipse, the centroid is located at 4∕3 with respect to the origin.
The static equilibrium equations yield a further unknown, which is the bending moment at the fixed end of the beam.As such, we introduce the differential equation of the elastic curve, such that where  is the elastic modulus of the material,  x x is the second moment of area of the cross section,  is the bending moment and  is the deflection of the beam.Equation 30 is integrated twice and solved for .A set of boundary conditions are defined, the maximum deflection of the beam at ∕ = 0, and the zero displacement point at  = 0.
The maximum deflection occurs midspan of the beam, whilst zero deflection occurs at the fixed ends.
The resulting integrals are solved analytically in Mathematica and numerically with Python.Results of the uniform loading case are verified versus commercial solvers SkyCiv and ClearCalcs.The procedure described here can be applied to different types of loading provided that their spanwise distribution is known.Further examples of loading can be found for example in Taylor and Hunsaker (2020) [45].
Lastly, the maximum bending stresses occur at the most distant point from the neutral axis of the beam cross section [50] and are given by Abel This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d where  max is the maximum bending moment and  max is the distance from the neutral axis to the outermost point of the beam.Because we use a symmetric NACA 0015, then  max = 0.075.
The coupling of the structural and hydrodynamic model is performed in Python.This allows for a structurally computationally inexpensive analysis and for a powerful evaluation tool that can be used to design large scale to full scale wave cycloidal energy converters.In the next section, we present the flow field around a laboratory scale cycloidal rotor to understand the effects of attached and separated flow conditions in the loading of the foil.We then assess the validity of the hydrodynamic model in these two flow regimes.In the remaining of the paper we use the coupled hydro-structural model that we introduced in §3 and §5 to determine the optimum rotor radius and span to achieve a balanced hydrodynamic and structural performance in wave design conditions and at different sea states.

Vorticity flow fields
To understand better the fluid dynamics and structural response of the single foil cyclorotor, we firstly study the flow field around a laboratory scale size rotor.The flow field is computed through two dimensional RANS simulations.This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352

P r e p r i n t n o t p e e r r e v i e w e d
The rotor radius is  = 0.3 m and the chord length of the foil is  = 0.3 m.The submergence of the rotor is  0 = −2.5.
The foil rotates clockwise and the wave direction is from left to right.The initial position of the rotor is defined as The attached flow case shows that the vorticity field is similar for the four azimuthal positions.Positive areas of vorticity emerge near the leading edge and are distributed along the inner side of the foil.In contrast, negative areas of vorticity are generated in the outer side of the foil, but these areas occur downstream of the mid chord.The wake of the foil follows the circumferential path of the rotor and it starts to dissipate approximately at a distance of 2.5 downstream of the foil.The circumferential path is denoted by a black dotted line in the figures.Noteworthy, the wake is fully dissipated at a distance of 4.This is approximately equivalent to a full cycle of rotation  and is the lifetime of the wake used in §3.
The vortical flow case shows a more dynamic flow field in the four azimuthal positions.At  = 0 • , a leading edge vortex (LEV) of positive vorticity and trailing edge vortex (TEV) of negative vorticity start to emerge both from leading and trailing edge, respectively.At  = 90 • , a large LEV is formed and convected along the suction side of the foil.The LEV forms a vortex pair with the coherent TEV that is located at a distance of about 1 downstream of the trailing edge.Subsequently, at  = 180 • , the LEV has dissipated and negative vorticity is shed in the form of a starting trailing edge vortex.A new layer of positive vorticity starts to form and to detach on the inner side of the foil.Lastly, at  = 270 • , the positive vorticity layer observed at  = 180 • has rolled up into a leading edge vortex, which induces a counter rotating secondary vortex of negative vorticity on the surface of the foil .At the same time, the TEV travels along the circumferential path and is found approximately at a distance of 2.5, and has started to dissipate.
It is noted that in the vortical flow condition, the LEVs convect along the suction side of the foil at about  = 0.5  and this will possibly be reflected in the instantaneous loading of the foil.Furthermore, the suction side is the inner side of the foil throughout the full cycle of rotation.This is similar to what is typically observed in a vertical axis wind turbine under dynamic stall, where the suction side is also located at the inner side of the foil [44,4].However, in VAWTs, the suction side remains on the inner side of the foil only for half a cycle of the rotation, whilst for the second This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d half the suction side changes to the outer side of the foil.For a cyclorotor in regular waves, the suction side remains on the inner side during the full rotation.This is because the wave velocity component acts always normal to the outer side of the foil.
In the next section, we select the instantaneous flow field at  = 90 • as a representative case of flow, for both attached and vortical conditions to analyse the flow field around the cyclorotor in more detail.

Topology of flow
The instantaneous streamlines surrounding the foils at  = 90 • for attached and vortical flow conditions are shown in figure 10a and figure 10b, respectively.This position is selected because it is when the foil is closest to the free surface, and the one where the foil experiences the highest angle of attack and therefore the highest loading.In the figures, the streamlines highlight the dominant directions of the flow around the foil.
Figure 10a shows that the dominant flow direction on the pressure side and upstream of the leading edge of foil is that of the wave velocity component.In these areas, the direction of the streamlines is mostly downwards.We recall that the wave velocity component acts normal to the motion of the foil because  = 90 • .In contrast, figure 10b shows that although the wave direction also points downwards, as evidenced behind the trailing edge of the foil at the upper left side of the figure, there is a strong jet of flow opposing the wave velocity component due to the vortex pair formed by the LEV and TEV.This opposite jet forces the flow on the pressure side of the foil to bend and be more tangential to the foil, as opposed to more normal, as it is in the case of figure 10a.
Two recirculation zones are identified in both of the figures with the streamlines.In the attached flow case of figure 10a, a small circulation bubble is present on the suction side of the foil.This recirculating zone is also common in steady translating flow on the concave side of curved surfaces, such as circular arcs at low Reynolds numbers [9].
A larger circulation zone with high content of positive vorticity is formed on the vortical flow case in figure 10b.In this case, an LEV with a diameter size of approximately 0.5 appears on the suction side of the foil.Contrary to the circulation zone in figure 10a, the LEV in figure 10b is mirrored by a TEV of opposite circulation.As such, according to impulse theory, the vortex pair will have an instantaneous influence in the force time stamp of the foil through growth rate of circulation ( Γ) and through the advection velocity of the LEV ( ḋ) [36,47].
Because the single point vortex model does not account for the flow physics that occur in the vortical flow case, it is likely that the analytical model might not reproduce the force signature under this scenario.This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d

Point vortex model assessment
In this section, we compare the forces computed with the point vortex model to those computed with RANS CFD simulations.The methodology for the CFD simulations is summarised in §4 and further developed in in Olbert et al. [35,34].The rotor parameters , ,  0 and  used in both CFD and point vortex model are the same as those described in section 6.1.The wave testing parameters   and   for each case are the same as those described in § 6.1.
We test the same two cases that we studied in the previous sections: attached and vortical flow cases.The average Reynolds number   for each of the two cases is computed by considering the average angular velocity (  ) of the foil.Here,   = (2∕  ).Hence,   = 300, 000 for the attached flow case, and   = 250, 000 for the vortical flow case.We note that although the Reynolds number is slightly different between the two cases, any significant difference in the loading is due to the different flow physics rather than by the effect of the Reynolds number.
In the single point vortex model, the range of  oscillations for the attached flow case ranges between 6 • <  ≤ 14 • , whilst for the vortical flow case, the range is 14   This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d to trough is greater in the vortical flow case than in attached flow case.The greater variance is due to the shedding of leading and trailing edge vortices in the wake of the vortical flow case, whilst the reduced variance of the attached flow case is, as shown in §6.1, due to the similarity of the flow throughout the full rotation of the rotor.
Although the point vortex model provides a satisfactory estimate of the mean value of the forces for the attached flow case, it overestimates the variance of the radial force for the attached flow case, whilst it underestimates the variance for the vortical flow case.The reason for the overestimation of the variance for the attached flow case is thought to be due to error in the estimation of the oscillations of the angle of attack and also, due to viscous effects around the foil that are not considered in the single point vortex model.In contrast, the underestimation of variance for the vortical flow case is due to the missing flow physics that the analytical model does not account for, i.e. leading and trailing edge vortices.Furthermore, the point vortex model does not account for dynamic stall effects, such as delay of separation and vortex lift [39,10].
However, because the point vortex model predicts the mean radial and tangential forces within reasonable accuracy with respect to the CFD computations, we consider the PV model to be a useful design tool, for attached flow cases, since it can predict the average mean forces that the foil experiences.
it is envisioned that in order for the rotor to balance hydrodynamic and structural performance, the rotor will need to operate Furthermore, the computational speed of the point vortex model (1 minute/case) compared to the time required to carry out one CFD simulation (24 hours/case) in a standard desktop computer, justifies the use of the point vortex model as a first design tool for cycloidal rotors.Furthermore, it allows for coupling of the structural model to allow an initial assessment of the structural design of this type of rotors.Lastly, because the rotor needs to be structurally resilient, it needs to operate in attached flow conditions.This is because in vortical flow conditions, similarly to vertical axis wind turbines [30], the variance of the forces can result in an increase in fatigue damage and therefore reduce the life of the rotor.As such, because the optimal operating conditions are those of attached flow, in §6.5 and §6.6, we perform the structural design of a large scale rotor using the single point vortex model to predict mean loads on the foil, and to size the rotor radius and span.

Instantaneous radial loading on foil
To further demonstrate that the optimal structural operation of the rotor is in attached flow conditions, in this section, we compare the instantaneous radial force  of the attached and vortical flow cases computed with CFD.We analyse  because it is the force component that is responsible for the bending moments and stresses on the foil, and because it is the force component that showed the most drastic difference in figure 11.
We plot the instantaneous  during one full cycle  after 35 in the polar plot of figure 12 This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d model for the case of attached flow.We note that  is defined at the quarter chord position of the foil and that it is positive inwards pointing towards the central shaft.In the figure the azimuthal angle denotes the rotor position and the radial coordinate denotes  in Newtons.
Figure 12 shows that the loading on the foil remains relatively constant for the case of attached flow.This is in agreement with the flow field observations presented in figure 9, in which the flow around the foil remains largely Although we only have information from four azimuthal positions of the flow field around the foil, at  = 0 • ,  = 90 • ,  = 180 • and  = 270 • , the trend in  can be explained through these flow field snapshots and the impulse of a vortex pair [8,36,47].As evidenced in figure 9, a clear vortex pair is formed at  = 90 • , whilst the pair is at its infant state at  = 0 • .From the flow field observations, it can be observed that between  = 0 • and  = 90 • , the circulation of the LEV and TEV grows, and that the convection velocity of the LEV is slower than the one of the TEV.
Hence, both terms, circulation growth and advection velocity of the LEV explain the increase in  of the foil [36,47].This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d Then, the LEV is likely to reach the advection velocity of the TEV and it dissipates in the wake.This probably explains the initial drop in .The flow is attached to the foil at  = 180 • , and because the initial vortex pair has dissipated, it is likely that the stable value of  at about 40 N is due to attached flow over the range of 115 • <  ≤ 225 • .Lastly, at  = 270 • , an LEV is present on the suction side of the foil, however, in this case the TEV has started to dissipate and is further downstream of the foil.Hence a vortex pair with cores of opposite and equal circulation is not formed.
However, the presence of the LEV influences and reduces any bound circulation on the foil [6,3,47].As such, the force decreases in the vicinity of this azimuthal angle.Lastly, the force recovers at  = 0 • , as the new vortex pair starts to emerge.
In summary, it can be seen in figure 12, that vortical flow is an undesirable condition in terms of load oscillations and structural loading.And therefore, a preferred operating point for the rotor is in attached flow conditions.Hence, in the next sections, we will size a large scale cyclorotor in terms of radius and span using the point vortex model.
We recall that the model provides mean radial and tangential mean loads that lie within 15% and 3% of the CFD simulations, respectively.Therefore, the potential flow model is a useful tool to design the rotor in attached flow conditions.Furthermore, because of the large oscillations that occur with vortical flow, it is expected that the optimal hydrodynamic and structural radius and span of the rotor will steer the rotor to operate in attached flow conditions and below the stall angle.This is similar to what has been observed in tidal turbines that operate in optimal conditions, in which attached flow dominates and dynamic stall is not influential [23,39].

Operational and structural design at wave design conditions
In this section we focus on finding the optimal operational phase and structural parameters for a large scale rotor under regular waves.Specifically, we focus on the selection of a radius that optimises power extraction and of a span length that maintains the bending stresses at the allowable stress level.
We first illustrate a sample case assuming wave design conditions, i.e.   = 9 s and = 1.2 m, as described in §2.
We plot the contour plots of the mean power output  and the maximum bending stresses  max in figures 13a and 13b, respectively.In both figures,  is plotted on the horizontal axis over a range of 0 • to 180 • , whilst  is plotted on the vertical axis over a range of 1 to 12 m.The foil has a chord length equal to .This is based on the findings of Siegel [40], who suggests optimal power capture when ∕ = 1.We find that for the wave design condition, a zero pitch angle on the foil yields negative power throughout most of the tested radii.As such, we apply a pitch angle of 5 • which is enough to increase the power output of the rotor but also, not enough to exceed the static stall angle.We note that this pitch angle is within the range of theoretical pitch angles that have been proposed for the operation of wave cycloidal rotors [41].This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d Figure 13a shows that for any , the optimal phase opt is  = 90 • .This is because at this , the tangential component of the lift force acting on the foils is maximised.Any other  results in a drop of the tangential force and therefore, of the mean output power.Noteworthy, when  = 0 • and  = 180 • , i.e. when the rotor is rotating in phase with the wave,  drops below zero.The figure also shows that  grows with  until an optimal  value, after which  starts to drop again.It can be seen that the optimal operating condition is when  is optimal and  = 90 • .For this particular case  opt is 8.5 m, as shown in figure 13a.
Figure 13b shows that  max grows with  and that the bending stresses remain mostly independent of  over a range of 60  Once  opt and  opt are determined, we can select the optimum span of the foil  opt .The larger the span, the more power the cyclorotor can produce, however, the loads and the span cannot be infinitely large due to their impact on the bending stresses of the foils.In this work, we apply the condition that the maximum bending stresses  max ≤   .We recall, from This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d  The span of the rotor is given by the intersection of the  opt curve (solid black line) and the   curve (dotted green line).At this intersection point, a vertical line is drawn to determine the span of the rotor that will satisfy ∕  .We find that for wave design conditions ( = 1.2 m and   = 9 s),  opt ∕ opt ≈ 3. Different materials for the foils could yield different ∕ ratios.For instance, the  opt ∕ opt ratio of the Atargis cyclorotor is  opt ∕ opt ≈ 10 with foils made of composite material [42].However, the methodology demonstrated here for sizing the rotor is independent of the material of the foil.In the next section, we explore the effect of different sea states in the  opt ∕ opt ratio.We note that although only results for uniform loading are presented, elliptical loading yields a reduction in  opt of about 5%.

Sizing of rotor for different sea states
The previous sections showed that there is an optimum radius  opt and optimum phase  opt at which the mean power output is maximised.We also showed that assuming  opt and  opt , the optimum span of the rotor  opt can be found, so that the maximum bending stresses at the fixed end of the foil remain at the allowable stress level   .In this section, we find the optimum radius  opt for different wave conditions assuming operation at  opt .We then find  opt for the same range of wave conditions.Results for  opt are shown in figure 15a. between 6 and 12 s, it is envisioned that for this particular site, the value of  opt would be between 5 to 10 m.This is agreement for example with the cycloidal wave rotor designed by Atargis corporation [42], which has been sized with a radius of 6 m.
Secondly, we find  opt for the same range of wave conditions tested in figure 15a.We note that we consider a chord length equal to  opt , as proposed by Siegel [41,42].Results for  opt are shown in figure 15b through the  opt ∕ opt ratio.
It can be seen that the shape of the  opt ∕ opt ratio mirrors the  opt matrix, and that  opt ∕ opt decreases from 4 to 1, as both   and   increase.Given that the probability of sea states is higher towards the bottom left side of figure 15b, we would expect a large scale rotor to to have an  opt ∕ opt ratio between 4 to 2. As an example,    of a 5 m rotor radius would be approximately 20 m.These dimensions can vary, for example, if different materials and allowable stress criteria than the ones used in this paper are considered.Here, we recall that we used offshore steel and   = 117 MPa, as indicated in table 5.In these results, we assume uniform loading, although results for elliptical loading yielded lower ratios by approximately 5%.
Finally, to assess the power capabilities of the single foil rotor and the stress level at the fixed end of the foils after sizing  opt and  opt , we compute the average power ( ) and the maximum bending stresses ( max ) in figure 16a and 16b, respectively.We use the same   and   combinations that we used in figures 15a and 15b.
We observe in figure 16a that  increases with   and   , whilst the stress level in figure 16b remains within ± 3% of 117 MPa.The increase in  with wave height is in agreement with previous power matrices of wave cyclorotors [42,18].We note however, that different to these studies, the power matrix presented here shows the maximum  at the top right corner of figure 16a and that power generation is possible also at very low   and   .This is because the power matrix of this work considers optimum power production through  opt , whilst at the same time, structural reliability by maintaining the stress level at the allowable threshold, as shown in figure 16b.
The results presented in this section highlight two important aspects in wave cyclorotor design that have not been addressed previously in the literature.Firstly, that through appropriate  and  sizing or that through variable  and , it is theoretically possible to generate power with wave cyclorotors even at low   and low   .This increases the range of locations where wave cycloidal rotors could be deployed.Secondly, that if a designer opts to increase the radius of the rotor, the span of the foil needs to be reduced to remain structurally resilient.However, this also means that a large span rotor could be reinforced with intermediate supports to reduce the free hanging parts of the foil.As such, our results provide guidelines for power production and resilient structural design at the significant wave conditions of a site.In addition to appropriate sizing, several control techniques, such as variable rotational velocity [17] and fatigue This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d damage mitigation strategies, such as passively pitching foils [7], could help in increasing power extraction, whilst reducing the loads on the foil and cyclorotor.

Conclusions
In this paper, the hydrodynamics and the structural design of a single foil wave cycloidal rotor in regular waves This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d due to fatigue damage.In contrast, vortical flows yield large radial amplitude oscillations, which are undesirable for the structural reliability of the rotor.
By comparing the CFD with the potential flow model results, we found that under attached flow conditions, the potential flow model yields reasonable accuracy, within 15% of the mean radial and within 3% of the mean tangential forces acting on the single foil.Contrarily, the analytical model underestimates the forces in the vortical flow case by at least 50%.This is likely due to the fact that dynamic stall is not accounted for.However, because the optimal hydrodynamic and structural operation of the rotor is expected to occur under attached flow conditions, the analytical model is considered useful and is utilised to to provide design guidelines by assuming optimal operating conditions for large scale rotors.
The design guidelines show that the optimal phase of operation is 90 • , and that there is an optimal radius at which the mean power production is maximised.The optimal radius can be used to determine the span of the foil that will keep the bending stresses at the allowable stress level.We find that because the optimal span to radius ratio changes This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d Abel Arredondo-Galeana et al.: Preprint submitted to Elsevier Page 2 of 32

Figure 1 :
Figure 1: LiftWEC wave bladed cyclorotor in operation near water surface and supported by two triangular frames.

Figure 2 :
Figure 2: Scatter plot showing energy period   versus significant wave height   of seasonal data of a point in the North Atlantic at the coast of France, located at 47.84 • N, 4.83 • W,

Figure 3 :
Figure 3: Side view of a single foil wave cycloidal rotor rotor showing the lift and drag forces (, ) on the hydrofoil, the wave velocity , the velocity due to the rotation of the hydrofoils , the relative velocity  and the operational phase  at ∕ = 0.

Figure 4 :
Figure 4: Angle of attack oscillations () in hydrofoil 1 for one period of revolution at design sea state conditions.

Figure 5 :
Figure 5: Volume fraction and grid resolution in CFD-simulation of single foil rotor in regular wave.

Figure 6 :
Figure 6: Domain setup employed in CFD simulations with inlet, outlet, top and bottom boundaries, and the damping zone.

Figure 7 :
Figure 7: a) Hydrofoil cross section of a curved NACA 0015 and b) simplified and straightened hydrofoil cross section composed of three hollow cross sections: triangular, rectangular and elliptical.

Figure 8 :
Figure8: Foil of cyclorotor subject to uniform and elliptical loading.The maximum distributed load is    and   for uniform and elliptical loading, respectively.The origin of the -axis is defined in the figure and  is the span of the foil.
Arredondo-Galeana et al.: Preprint submitted to Elsevier Page 15 of 32

3 .
figure 3. Two cases are studied.One where the flow remains attached (  = 0.253 m,   = 1.829 s) to the foil and one where the flow separates and is dominated by vortex flow (  = 0.31 m,   = 2.462 s).These two cases are chosen because of the contrasting flow features that are likely to yield different structural loads.We refer to the cases as attached and vortical flow cases, respectively.Figures 9a-d and figures 9e-h show the vorticity flow fields for the attached and vortical flow conditions, respectively.Four azimuthal positions are shown:  = 0 • , 90 • , 180 • and 270 • , which are shown from left to right in the two rows of the figure 9.

Figure 10 :
Figure 10: Velocity fields around hydrofoil at  = 90 • for a) attached and b) vortical flow conditions.
estimation by the potential flow model of about 2 • .The resulting forces computed with CFD and with the point vortex model (PV) on the single foil of the cyclorotor are shown in figure11aand figure11bfor the attached and vortical flow cases, respectively.We use whisker plots to show the force data.The whisker plots show the mean value, the maximum and the minimum values of the radial and tangential forces during one cycle of rotation  after 35 .This ensures that any transient force resulting from initial accelerations on the CFD simulations are not present in the force results.

Figure 11a shows thatFigure 11 :
Figure11ashows that in the attached flow case, the mean values of the CFD simulations and the single point vortex model lies within 15 %, whilst the agreement of the mean tangential forces lies within 3 %.In contrast the agreement of the mean radial and tangential forces decreases in the vortical flow case.In this case, the mean value of the radial forces is underpredicted by the point vortex model by approximately 50%, whilst the mean value of the tangential forces is under predicted by approximately 80%.The CFD results also show that the oscillation of the radial force from peak . Results are shown for both attached and vortical flow cases.For comparison, we also present the average  computed with the point vortex Abel Arredondo-Galeana et al.: Preprint submitted to Elsevier Page 21 of 32 unaltered at  = 0 • , 90 • , 180 • and 270 • .The figure confirms that the average value of  predicted by the single point vortex model ( -blue dotted line) lies within 15% of the magnitude of  computed with CFD (black line).In contrast, for the vortical flow case,  (red line) increases from  = 0 • to  = 68 • to a maximum value of about 80 N.Then,  remains relatively constant until about  = 80 • , after which  drops and stabilises at around 115 • at approximately 40 N.Then,  remains stable at a value close to 40 N until 225 • , after which it drops further to its lowest value of 18 N at approximately 300 • .Lastly,  slowly recovers to a value of 40 N at 0 • to start the cycle again.

Figure 12 :
Figure 12: Radial forces  for attached and vortical flow conditions computed from CFD and from analytical model for the case of attached flow, during one period of revolution ∕  .

Figure 13 :
Figure 13: a) Mean output power  and b) maximum bending stresses  max plotted as a function of  and  for the design sea state conditions (  = 1.2 m and   = 9 s).

Figure 14 shows
Figure 14 shows  max versus different foil spans.We present results for uniform loading.In the plot,  max is plotted for a range of different radii ranging between 1 m to 12 m.The black line shows the maximum bending stresses in the foil at  opt , whilst the blue and pink line show the stresses at  min = 1 m and  max = 12, respectively.Stresses for intermediate radii are plotted with dashed gray lines in steps of 1 m.The allowable stress level   and the yield stress level   are shown with horizontal green and red dotted lines, respectively.

Figure 14 :
Figure 14: Maximum bending stresses ( max ) versus span () for different rotor radii.The blue, black and pink lines correspond to  min ,  opt and  max .The dotted gray lines correspond to intermediate radii.The green and red horizontal lines correspond to   and   , respectively.

Figure 15 :
Figure 15: a) Optimum radius  opt and b)  opt ∕ opt ratio for different wave conditions for a single foil cyclorotor have been studied.We do this by developing a potential flow model coupled to a structural model to design single foil cyclorotors in an efficient manner.The potential flow model considers unsteady wake effects and includes a correction for flow separation.To understand better the flow physics of the wave cycloidal rotor and the limitations of the potential flow model, we analysed the flow field around a single foil rotor with the aid of two dimensional RANS simulations.A laboratory scale type of rotor was studied.Two flow conditions were analysed, attached and vortical flow conditions.It was found that attached flow conditions are desirable to minimise radial loading fluctuations, and therefore reducing the potential Abel Arredondo-Galeana et al.: Preprint submitted to Elsevier Page 27 of 32

Figure 16 :
Figure 16: a) Power matrix in kW and b) stress level at fixed end of the foils with uniform loading in MPa for rotor at  opt ,  opt and  opt ∕ opt .

Abel
Arredondo-Galeana et al.: Preprint submitted to Elsevier Page 28 of 32This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d with different wave conditions, the dimensions for a fixed span fixed radius rotor are likely to be determined by the sea state with highest probability of occurrence.The novelty of the results presented in this paper include, for the first time in the literature of wave cycloidal rotors, a detailed flow field characterisation of a single foil rotor.We show, for the first time, under what flow conditions a single point vortex potential flow model can be used to estimate the forces of the rotor subject to regular waves.These findings contribute and pave the way to further advance the research and development of cycloidal rotors for wave energy conversion.
where  0 is the submergence of the rotor, () is the angular position measured with respect to the negative horizontal axis and positive clockwise.The horizontal and vertical velocity components of the point vortex (  ,   ) are   =  sin(()) Abel Arredondo-Galeana et al.: Preprint submitted to Elsevier Page 6 of 32This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn. 12) 0 = 2 and  0 = 0 • .Both values are based on the lift curve slope of a NACA 0015 foil at  ≥ 250, 000.The latter covers the range of  of the test cases presented subsequently in §6.Equation 19 is solved for  .In the equation,   is computed with static   and   data from a NACA 0015 at the  at which each test of this paper is studied.
Abel Arredondo-Galeana et al.: Preprint submitted to Elsevier Page 9 of 32This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d

Table 1 :
Mechanical properties of moderate strengthoffshore steel table, the allowable stress level is defined as one third of the yield stress level (  ).However, different thresholds can be selected according to design requirements.
In the next section, we investigate to what extent is the point vortex model accurate to predict loading in the case of attached and vortical flow scenarios.
Abel Arredondo-Galeana et al.: Preprint submitted to Elsevier Page 18 of 32 and Klimas, the stall angle for this foil at this  is   = 12 • .Because in the attached flow case, the computed  oscillations with the potential flow model exceed   by 2 • and stall does not occur, this could mean an error in the • <  ≤ 21 • .Because   and   data for a curved foil is not available in the literature, we consider stall angle   of a NACA 0015 at  = 360, 000 as a reference point.According to Sheldahl Abel Arredondo-Galeana et al.: Preprint submitted to Elsevier Page 19 of 32This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d • ≤  ≤ 120 • .The growth of bending stresses with  is due to the fact that  grows with .As such, event if the static stall angle is exceeded ( ≥  opt ) and  and  drop,  continues to increase due to an increase in .Hence, we utilise figure13ato determine  opt and  opt , and 13b to ensure that  max ≤   .
Abel Arredondo-Galeana et al.: Preprint submitted to Elsevier Page 24 of 32 Even though  opt changes with   and   , as shown in figure15a, it is expected that a large scale rotor will have a fixed radius.As such,  opt is likely to be determined by the sea state condition with the highest probability of occurrence in a given location.Recalling that in figure2, the sea state with highest occurrence are those with   ≤ 2 m and The horizontal axis shows   over a range of 6 to 16 s and the vertical axis   over a range of 1 to 5 m.The figure shows that  opt grows with   and   .Abel Arredondo-Galeana et al.: Preprint submitted to Elsevier Page 25 of 32This preprint research paper has not been peer reviewed.Electronic copy available at: https://ssrn.com/abstract=4239352P r e p r i n t n o t p e e r r e v i e w e d