Control of an oscillating-water-column wave power plant for maximum energy production

Applied Ocean Research 24 (2002) 73–82 www.elsevier.com/locate/apor

Control of an oscillating-water-column wave power plant for maximum energy production A.F. de O. Falca˜o Department of Mechanical Engineering, Instituto Superior Te´cnico, Avenida Rovisco Pais, 1049-001 Lisbon, Portugal Received 3 May 2002; accepted 8 July 2002

Abstract A stochastic model was applied to devise an optimal algorithm for the rotational speed control of an oscillating-water-column (OWC) wave power plant equipped with a Wells turbine and to evaluate the average power output of the plant. The hydrodynamic coefficients for the OWC are assumed known (as functions of frequency), as well as the turbine performance curves. The whole model is based on linear control theory of a stochastic process, it being assumed that the sea surface elevation has a Gaussian probability density function. The optimal control law is expressed in terms of a simple relationship between the instantaneous values of the electromagnetic torque (to be applied on the generator rotor) and the rotational speed. It is remarkable that the optimal control algorithm was found to be practically insensitive to wave climate. A simple additional algorithm, accounting for constraints imposed by the electrical grid on power oscillations, was derived in order to complement the optimal control law. q 2002 Elsevier Science Ltd. All rights reserved. Keywords: Wave energy; Oscillating-water-column; Wells turbine; Control; Stochastic methods; Probabilistic methods; Optimization

1. Introduction The oscillating-water-column (OWC) is probably the most highly developed type of wave energy converter and one of the very few to have reached the stage of full-sized prototype [1,2]. The present paper addresses the control of OWC plants, especially on what concerns the rotational speed of their air turbines. In terms of time-averaged values (denoted by an over bar), the electrical power output from an OWC device may
w 2 L;
where Ww is the power absorbed be written as P
¼ W from the waves and L are the losses that occur in the energy conversion chain. The most significant losses are: (i) L1 due to real (viscous) fluid effects in the water; (ii) aerodynamic losses L2 in the turbine (including connecting ducts) and valves; (iii) bearing friction losses L3; (iv) electrical losses L4. From the energy production point of view, the plant should be controlled to maximize the value of P
in every sea state. This is to be done by controlling the turbine rotational speed (or alternately the torque applied by the electrical generator upon the turbine) and also possibly by controlling the valve system. Phase control by latching or by variableE-mail address: [email protected] (A.F. de O. Falca˜o).

geometry turbine (e.g. variable pitch Wells turbine) is not considered here. It is known that, on the one hand, the turbine aerodynamic efficiency depends on flow-rate and rotational speed. On the other hand, variations in the turbine rotational speed affect the flow-rate versus pressure-head curve (i.e. alter the turbine damping), and, in this way, modify the hydrodynamic radiation admittance [3,4]; this in turn affects the amount of energy absorbed from the waves. For these reasons, the hydrodynamics of the wave energy absorption and the turbine aerodynamic performance are coupled through the turbine rotational speed and should be considered jointly when optimizing the plant control algorithm. The problem of the rotational speed control of an OWC in random waves was dealt with by Justino and Falca˜o [3], who considered, numerically simulated and compared several control algorithms, including their stability. The control strategy proposed in Ref. [3] as the most suitable one is based on dimensional analysis applied to turbomachinery. If the flow is assumed incompressible and the effects due to Reynolds number and Mach number variation neglected, the instantaneous value of the aerodynamic efficiency of a given turbine is known to be a function only of Tt =N 2 ; where Tt is the aerodynamically produced torque on the turbine rotor

0141-1187/02/$ - see front matter q 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 1 4 1 - 1 1 8 7 ( 0 2 ) 0 0 0 2 1 - 4

74

A.F. de O. Falca˜o / Applied Ocean Research 24 (2002) 73–82

and N is rotational speed. If the bearing friction torque is ignored, the time-averaged values of Tt and Te (electromagnetic torque on the generator rotor) must be equal. This was used in Ref. [3] to propose, as control algorithm, a relationship of the type Te ðNÞ ¼ constant £ N 2 to be implemented in the programmable logic controller (PLC) of the plant. This approach aimed at keeping the timeaveraged value of the turbine aerodynamic efficiency at its maximum, but ignored the effect of varying rotational speed upon the hydrodynamics of wave energy absorption, and did not propose any non-empirical method for determining the value of the proportionality constant. In the present paper we look for a control algorithm to be implemented in the plant’s PLC relating the instantaneous electromagnetic torque (applied upon the generator rotor) to the instantaneous rotational speed; the algorithm should take into account the turbine efficiency as well as the hydrodynamics of wave energy absorption. In the optimization only losses L2 and L3 are considered. Viscous losses in the water, L1, are likely to be significant, but are difficult to model theoretically and are not expected to depend strongly on plant control. Electrical losses L4 are relatively small and their variation with generator rotational speed may reasonably be neglected. The existence of a valve is considered, its role being to control the air flow rate through the turbine in order to avoid aerodynamic stalling at the rotor blades, but not phase control (by a valve system or otherwise). The theory is based on a stochastic modelling of the OWC performance as developed in Ref. [4]. The statistical theory of random surface waves is briefly presented in Section 2. It is assumed that the local wave climate may adequately be represented by a set of sea states of known power spectrum and frequency of occurrence, the surface elevation being a Gaussian random variable for each sea state. The air pressure inside the chamber is regarded as the response of the assumedly linear system (consisting of the structure, chamber and Wells turbine) to the random waves. In particular, it is shown how the variance and the probability distribution of the chamber air pressure can be computed, it being supposed that the hydrodynamic coefficients of the device are known as functions of frequency (and possibly wave direction), and that the performance curves of the Wells turbine are given. A proportionality relationship (at constant rotational speed) between turbine flow-rate and pressure-head is a basic assumption underlying the whole analysis. Average values for the power output can then be easily obtained from the turbine instantaneous performance curves and the pressure probability density function. Dimensionless values are widely used for the turbine variables in order to extend the applicability of the results. It is shown how the stochastic method can be employed to optimize the rotational speed for maximum energy production. In Section 3, the theory is applied to the rotational speed control of the Wells turbine that equips the OWC plant built on the Island of Pico, Azores.

2. Theory We assume that the local wave climate may be represented by a set of sea states, each being a stationary stochastic ergodic process. For a given sea state, the probability density function f(z ) of the surface elevation z at a fixed observation point is supposed to be Gaussian, an assumption widely accepted in ocean engineering [5 – 7] and adopted in Ref. [4], and so we may write ! 1 z2 ð1Þ f ðzÞ ¼ pffiffiffiffi exp 2 2 ; 2sz 2psz where s2z is the variance, and sz is the standard deviation, of z. We assume here one-directional incident waves (multidirectional waves are dealt with in Ref. [4]). The variance is related to the (one-dimensional) spectral density by [4] ð1 s2z ¼ Sz ðvÞdv; ð2Þ 21

where v is radian frequency. We note that Sz ðvÞ is a twosided spectral density (per unit radian frequency bandwidth) defined in the interval ð21; 1Þ and is an even function. It is related to the more frequently used one-sided spectral power ð1Þ ð1Þ density Sð1Þ z ðvÞ by Sz ðvÞ ¼ 2Sz ðvÞ if v . 0; Sz ðvÞ ¼ 0 if v # 0: It is convenient to introduce dimensionless variables to characterize the turbine aerodynamic performance, and so we write

C¼

p ; ra N 2 D 2

ð3aÞ

F¼

m _ ; ra ND3

ð3bÞ

P¼

Pt : ra N 3 D5

ð3cÞ

Here C is non-dimensional pressure-head, F is nondimensional flow-rate, P is non-dimensional turbine power output, ra is atmospheric air density, m _ is mass flow-rate through the turbine, p is air pressure oscillation in the chamber, N is rotational speed (radians per unit time), Pt is turbine power output (bearing friction losses are ignored) and D is turbine rotor diameter. Note that here the losses in the connecting ducts are considered included in the turbine aerodynamic losses. If the effects of the variations in Reynolds number and Mach number are ignored, dimensional analysis allows us to write [4,8]

P ¼ fP ðCÞ;

ð4aÞ

F ¼ fQ ðCÞ;

ð4bÞ

the functions fP and fQ depending on the shape of the turbine but not on its size or rotational speed. It was found experimentally (with some theoretical foundation) for the Wells turbine that Eq. (4b) may be replaced approximately

A.F. de O. Falca˜o / Applied Ocean Research 24 (2002) 73–82

by

F ¼ K C;

ð5Þ

where K is constant for a given turbine geometry (is independent of turbine size or rotational speed); the adequacy of this assumption was discussed in detail in Ref. [4]. In what follows we take fP to be an even function (symmetric turbine). If linear water wave theory and Eq. (5) are adopted, and in addition the spring-like effect of air inside the chamber is modelled by a linearized isentropic relationship between air pressure and density, then it may be found [3,4,9] that the system, whose input is incident wave elevation and output is inside air pressure, is a linear one. It follows that the chamber air pressure oscillation p(t ) has a Gaussian distribution [4,10], and its variance s2p is related to Sz ðvÞ by ð1 s2p ¼ 2 Sz ðvÞlGðvÞLðvÞl2 dv: ð6Þ

75

same applies to Pe, which allows us to write Pe ¼ P
t 2 L
3 ; where L
3 is the averaged value of the bearing friction loss. Let us consider a given sea state and ignore, for the moment, the bearing friction loss. From the definitions (10a) and (10b), we easily obtain 1 dP
t 3P
¼ þ 3 sC D sp dN

! N dsp dP
22 : sp dN d sC

ð11Þ

The condition for maximum turbine power output implies dP
t =dN ¼ 0; and so, from Eq. (11), it follows that

sC dP
¼ P
dsC

3 : N dsp 22 sp dN

ð12Þ

0

Here G(v ) is the excitation-volume-flow coefficient (as defined in Ref. [9]), and    21 KD v V0 þB þi L¼ þC ; ð7Þ ra N gpa where pa is the atmospheric pressure, g ¼ cp =cv is the specific heat ratio for air, V0 is the undisturbed value of the air chamber volume, and the hydrodynamic coefficients B(v ) and C(v ) are the radiation conductance and the radiation susceptance, respectively [9]. The averaged value of the turbine power output (bearing friction losses ignored) is !   3 5 ð1 2 r N D p p a P
t ¼ pffiffiffiffi exp 2 2 fP dp; ð8Þ 2 sp ra N 2 D2 2psp 21 or in dimensionless form 1 P
¼ pffiffiffiffi 2psC

ð1

2

P
t ¼ ra D5 P
p N 3 :

!

C exp 2 2 fP ðCÞdC; 2 sC 21

P
t ; ra N 3 D 5 sp sC ¼ : ra N 2 D 2

ð13Þ

ð9Þ

where

P
¼

For a given turbine, the left-hand-side of Eq. (12) can be regarded as a known function of sC which we denote by uðsC Þ: We recall that, for a given sea state (characterized by its spectral distribution), Eqs. (6) and (7) give sp, and hence the right-hand-side of Eq. (12), as functions of N. It now becomes clear that Eq. (12) yields the value of N that maximizes the averaged turbine power output P
t : The optimal control strategy to be implemented in the PLC of the plant consists essentially in setting Pe ¼ P
t ; P
t being related to N through Eq. (12). If the effect of varying rotational speed upon the wave energy absorption is neglected, as was done in Ref. [3], then dsp =dN ¼ 0 and Eq. (12) reduces to uðsC Þ ¼ 3=2: We denote the solution of this equation by spC and the corresponding value of P
by P
p : From Eq. (10a) we find

ð10aÞ ð10bÞ

We note that the function fP is assumed known for the turbine under consideration (from testing at reduced or full scale, or from numerical modelling), and, from Eq. (9), the same can be said of the function P
ðsC Þ: The problem to be solved is to find the optimal relationship, for each sea state, between the electromagnetic torque Te (or the corresponding power Pe ¼ Te N) and the rotational speed N. The moment of inertia of the rotating parts is supposed large, so that the variations in N are relatively small in a given sea state. We assume that the

We note that P
p is a characteristic of the turbine that does not depend on rotational speed or on sea state. If the electric power Pe is equated to P
t ; Eq. (13) reproduces the control strategy Pe ðNÞ ¼ constant £ N 3 proposed in Ref. [3] and, in addition, supplies the value of the proportionality constant. If the effect of varying rotational speed upon the wave energy absorption is to be considered, then the term ðN=sp Þ dsp =dN cannot be neglected in Eq. (12). Since the air pressure standard deviation sp depends on N and also on the spectral distribution that characterizes the sea state under consideration, the relationship (to be implemented in the plant’s PLC) between Pe and N is expected to become more complicated. However, it will be shown in Section 3 (where flow control by a relief valve is considered and bearing friction losses are accounted for) that this may be replaced with good approximation by a simple relationship between average power output P
t and rotational speed N, regardless of the wave climate.

76

A.F. de O. Falca˜o / Applied Ocean Research 24 (2002) 73–82

Fig. 1. Local wave climate, represented by 44 sea states (Te, Hs). The area of each circle is proportional to the contribution of the sea state to the total annual wave energy.

3. Example calculation: the Pico plant 3.1. The plant The method developed above is now applied to the control of the OWC shoreline wave power plant constructed on the Island of Pico, Azores, Portugal. The plant, described in Ref. [1], has a concrete structure (of 12 m £ 12 m internal square cross-section at free-surface level) that spans a natural gully whose depth is about 8 m. It is equipped with a horizontal axis Wells turbine driving a variable-speed electrical generator. Results from wave measurements at the plant’s location, in conjunction with longer-term results from forecast/hindcast numerical modelling for the same ocean area, were used to establish a set of 44 sea states characterizing the local wave climate. Each sea state is represented by its significant wave height Hs, its energy period Te and its frequency of occurrence f (see the table in Ref. [11]). Fig. 1 shows a plot of Hs versus Te for the 44 points, with a circle

Fig. 2. Hydrodynamic coefficients for the Pico plant: radiation conductance B (solid line), radiation susceptance C (dashed line) and excitation-volumeflow coefficient G (chain line).

Fig. 3. Dimensionless plot of Wells turbine power output versus pressurehead: P(C ) (solid line), P
ðsC Þ without valve control (dashed line), and P
ðsC Þ with valve-controlled flow (chain line).

whose area is proportional to the sea state contribution to the total annual energy flux centred at each point. As in Ref. [4], we adopt the following standard formula (Ref. [7], p. 28) for the frequency spectrum of each sea state Sz ðvÞ ¼ 131:5Hs2 Te24 v25 expð21054Te24 v24 Þ;

ð14Þ

where Hs is the significant wave height and Te is the energy period. Curves for the hydrodynamic coefficients lGðvÞl; B(v ) and C(v ) are given in Fig. 2. They were obtained from numerical values computed for the Pico plant by Brito-Melo et al. [11] with a three-dimensional boundary element code. The complex shape of the curves reflects the irregular bottom and surrounding walls modelled. The Wells turbine has an 8-bladed rotor of D ¼ 2.3 m outer diameter and Di ¼ 1:36 m inner diameter, with a row of guide vanes on each side of the rotor. The chord of the rotor blades is constant c ¼ 0:375 m and its cross-section varies continuously from NACA 0015 at the hub to NACA 0012 at the tip. The curve P ¼ fP ðCÞ can be found in Ref. [3] and is represented by the solid line in Fig. 3. The curve F ¼ fQ ðCÞ was approximated by the straight line F ¼ K C; with K ¼ 0:6803 [3]. It may be seen that P increases with C, and reaches a maximum Pmax ¼ 0:00213 at a critical value Ccr ¼ 0:067; above which P dramatically drops due to aerodynamic stalling at the rotor blades. Eq. (9) may be used to compute the average power output P
in random waves as a function of the r.m.s. value, sC, of the dimensionless pressure-head C. In the computation fP ðCÞ was taken as an even function. Since the integration in Eq. (9) extends to infinity, the curve P(C ) in Fig. 3 was extended by setting P ¼ constant ¼ 0:00074 for C $ 0:095; which is probably a conservative estimate of the turbine performance under stalled conditions. The dashed curve in Fig. 3 represents P
versus sC. The solid and the dashed curves are very different from each other, which is not surprising since the former represents instantaneous, and the latter averaged, values of (dimensionless) power output. Besides, the curve of P
ðsC Þ (unlike P(C )) shows relatively little variation to the right of its maximum point.

A.F. de O. Falca˜o / Applied Ocean Research 24 (2002) 73–82

Fig. 4. Time-averaged bearing friction power loss N T
a ; due to axial load, versus sC, for N ¼ 157.1 rad s21 (1500 rpm), without (solid line) and with relief valve (dashed line).

3.2. Control valves A way of preventing the turbine rotor blade stalling (and avoiding the drop in power output that occurs whenever lCl exceeds its critical value Ccr) is to control the air flow rate through the turbine by installing either a valve in series or a valve in parallel with the turbine. This has been examined in Ref. [12] where a detailed comparative analysis has been carried out in the time domain. It should be pointed out that partially closing a control valve in series or, alternately, opening a control valve in parallel, introduces a nonlinearity in the relationship between the air flow rate leaving (or entering) the chamber and the excess air pressure p in the chamber. Consequently the whole system ceases to be linear and the probability density function of the air pressure in the chamber is no longer Gaussian. In order to keep the simplicity of the present approach, we have assumed that the introduction of a properly controlled valve system is equivalent, for the purpose of computing the plant performance, to replacing the turbine power curve, P ¼ fP ðCÞ; for lCl $ Ccr ; by lPl ¼ constant ¼ Pmax ; this approximation is discussed in detail in Ref. [4]. The chain line in Fig. 3 represents the corresponding results for the average power P
ðsC Þ; and shows that P
increases with sC, tending asymptotically to Pmax as should be expected. (We note that a linear relationship, F ¼ K C; is still assumed between the dimensionless excess pressure in the chamber C and the dimensionless volume flow rate F through the turbine.)

77

on the atmospheric air pressure pa (Eq. (7)). We will assume that pa is approximately constant and equal to the standard atmospheric pressure at sea level pa ¼ 1:013 £ 105 Pa: Then, for a given sea state (characterized by its spectral density Sz) and given chamber geometry, sp may be regarded as a function of KX that we denote by sp ¼ Fp ðKXÞ ¼ Fp ðKDðra NÞ21 Þ: Since the dimensionless timeaveraged power output of the turbine P
is a known function of sC (Fig. 3), it can be regarded (for fixed ra and D ) as a function of N. The time-averaged turbine power can be calculated from P
t ¼ ra N 3 D5 P
: We note that Pt ¼ NLt ; where Lt is the aerodynamic torque on the turbine rotor. For constant rotational speed, the turbine net-power is P
n ¼ P
t 2 L
3 ; where L3 is the (turbine and generator) bearing friction loss. 3.4. Bearing friction losses We recall that the rotational axis is horizontal, and note that the loads on the bearings are: (i) the weight of rotating parts; (ii) an aerodynamically produced axial force on the turbine rotor. The contribution Tr to the bearing friction torque due to the radial load (weight of the rotating parts) could be calculated from the record of the rotational speed decay with time at zero electrical power and zero air flow rate (below about 100 rpm the aerodynamic drag becomes negligible), and was found to be approximately constant Tr ¼ 12:8 N m: The aerodynamically produced axial force depends on the instantaneous flow-rate through the turbine; this is difficult to measure in practice, but can be approximately computed from theoretical flow modelling. The expressions derived in Appendix A may be used to evaluate approximately the axial load Fa on the bearings (in the present case it is n ¼ D=Di ¼ 0:591 and a1 ¼ 62:58). In order to compute the corresponding friction torque Ta, we assume that Ta =lFa l ¼ Tr =Mg ¼ l; where M ¼ 2850 kg is the total mass of the rotating elements, and obtain l ¼ 0:458 £ 1023 m: We take the resulting friction torque to be simply Tr þ Ta ; and write for the lost power due to bearing friction L3 ¼ NðTr þ Ta Þ: Fig. 4 shows a plot of N T
a (average power loss due to axial load) versus sC for N ¼ 157.1 rad s21 (1500 rpm). Especially for the higher values of sC, the pressure difference across the turbine (and hence the axial load) is reduced by the action of the relief valve (if there is one), and the same is true for N T
a as shown in Fig. 4 (compare the solid and broken lines). 3.5. Numerical results

3.3. Plant performance in random waves The variance sp of the air pressure in the chamber may be calculated from the incident wave spectral density by using Eqs. (6) and (7). We note that the excitation-volumeflow coefficient G depends only on the chamber geometry, whereas L depends also on KX (where X ¼ Dðra NÞ21 ) and

In the computations of the turbine power output, the following values were taken: g ¼ 9.8 m s22, rw ¼ 1025 kg m23 ; V0 ¼ 1050 m3 ; g ¼ 1:4; ra ¼ 1:25 kg m23 : Since we are interested here in the rotational speed control, we computed, for each sea state, the optimal rotational speed, i.e. the value of N that maximizes the

78

A.F. de O. Falca˜o / Applied Ocean Research 24 (2002) 73–82

Fig. 5. Time-averaged net-power from turbine versus rotational speed (log – log plot), under optimal conditions, without flow control by valve. Each point is represented by a circle whose area is proportional to the contribution of the corresponding sea state to the total annual wave energy. The curve’s equation is P
n ¼ 1:583 £ 1025 N 3:1616 ; with P
n in kW and N in rad s21.

Fig. 7. As in Fig. 5, for a plant with flow control by relief valve. The curve’s equation is Pe ¼ P
n ¼ 1:093 £ 1025 N 3:2361 :

time-averaged value of the turbine net-power output P
n ¼ P
t 2 L
3 : Fig. 5 shows a log – log plot of P
n versus N under optimal conditions, without any flow control by relief valve. As in Fig. 1, the area of each circle is proportional to the contribution of the corresponding sea state to the total annual wave energy. It is striking that, unlike in Fig. 1 (where the points are widely scattered), the points, in Fig. 5, are almost perfectly aligned along a single curve; its equation was found to be P
n ¼ 1:224 £ 1025 N 3:1563 (P
n in kW, N in rad s21). It is appropriate at this stage to discuss how this curve can be used as a control law. We begin by noting that the set of 44 sea states was taken as a convenient representation of the local wave climate. It could as well be replaced by a much larger set of much shorter ‘sea states’ each with a duration of a few wave periods along which the inertia of the rotating parts would ensure the rotational speed not to vary significantly. It is to be expected that the points representing the optimal plant performance under this larger set of

shorter sea states would still be aligned along the same curve in an ðN; P
n Þ plot. If this optimal relationship between average net-power output P
n and rotational speed N is to be kept, then one should set Pe ð¼ NLe Þ ¼ P
n ; where Le is the electromagnetic torque applied on the generator rotor and Pe the corresponding power. In the case under consideration here, the optimal control law should be Pe ¼ 1:224 £ 1025 N 3:1563 : The values of N for the points representing the more energetic sea states in Fig. 5 considerably exceed the maximum allowable value for the Pico turbine, that has been specified by the manufacturers as about 1500 rpm (157.1 rad s21) (the constraints concern centrifugal stresses at the turbine rotor, Mach number effects, and electrical equipment performance). Fig. 6 (linear plot) is a modified, and more realistic, version of Fig. 5, in which the rotational speed is required to satisfy N # 157:1 rad s21 : It may be seen that the introduced constraint substantially reduces the turbine power in the more energetic sea states, as should be expected (the highest plotted value of P
n is now 188 kW as compared with 1135 kW). The same kind of results are shown in Figs. 7 and 8 for

Fig. 6. As in Fig. 5 (linear plot), but with rotational speed constraint N # 157:1 rad s21 :

Fig. 8. As in Fig. 7 (linear plot), but with rotational speed constraint N # 157:1 rad s21 :

A.F. de O. Falca˜o / Applied Ocean Research 24 (2002) 73–82

79

of a properly controlled relief valve may provide, in the case under consideration, an increase in annual energy production of about 37%. It may be of interest to compare the turbine output if the cube-law were applied instead of the optimal control algorithm. We assume a relief valve to be used, and consider only the 14 (less energetic) sea states for which the rotational speed is not constrained by N # 157:1 rad s21 (they represent 62.7% of the total time). We find, in average, 58.8 kW instead of 62.7 kW, i.e. a reduction of 5.2%. Correspondingly, the rotational speed is smaller (as should be expected from the comparison of the curves in Fig. 9), with ratios between 0.81, for the less energetic sea state, and 0.89, for the most energetic one. Fig. 9. Comparison of rotational speed controls laws. Thick (thin) lines concern situations without (with) control valve: Pe ¼ 1:224 £ 1025 N 3:1563 ( ¼ 1.093 £ 1025N 3.2361). Dashed lines represent simplified laws Ppe ¼ 3:344 £ 1025 N 3 ( ¼ 5.054 £ 1025N 3).

the case when a controlled relief valve prevents the air flow rate through the turbine from exceeding the stall-free limit. The new control equation is now Pe ¼ P
n ¼ 1:093 £ 1025 N 3:2361 : Comparing Figs. 6 and 8 (for which it is N # 157:1 rad s21 ), it may be seen that the use of a relief valve may increase very substantially the power delivered by the turbine in the more energetic sea states: the highest plotted value of P
n in Fig. 8 is 451 kW, and only 188 kW in Fig. 6. It may be of interest to compare the optimal control curve Pe ¼ P
n ¼ 1:224 £ 1025 N 3:1563 ; shown in Fig. 5, with the results from the simplified control law given by Eq. (13). We find spC ¼ 0:03253 and P
p ¼ 0:0004157; from which we obtain Ppe ðNÞ ¼ 3:344 £ 1025 N 3 : If a relief valve is operating, the corresponding control laws are Pe ¼ P
n ¼ 1:093 £ 1025 N 3:2361 and Ppe ðNÞ ¼ 5:054 £ 1025 N 3 : Results from these four equations are plotted in Fig. 9. The thick lines correspond to the situation without relief valve whereas the thin lines mean that the flow is controlled by a relief valve. The dashed (thick and thin) lines represent the simplified control laws. It may be seen that the N 3-laws (dashed lines as compared with solid lines) specify significantly larger electrical power (and electromagnetic torque) for the same N, although (since they are not optimal laws) it is obvious that they result in less produced energy. The annual-averaged net-power output from the turbine may be computed from P
n;annual ¼

44 X

fj P
n;j :

j¼1

Its value was found to be equal to 100.3 kW if no relief valve is used, and 137.0 kW if the air flow rate is prevented from exceeding the stall-free limit by a relief valve, in both cases with the rotational velocity restricted to a maximum value of 157.1 rad s21 (1500 rpm). This means that the use

3.6. Constraints introduced by the electrical grid The oscillations in electrical power output from an OWC wave energy plant can easily be absorbed by a large electrical grid, but may introduce unacceptable disturbances into a small isolated grid, as is typically the case of Pico Island. The capability of storing kinetic energy in, and releasing it from, a flywheel is a way of smoothing the oscillations in electrical power delivered to the grid. However, this has to be associated with a control law that allows the rotational speed to oscillate. For this reason, one should avoid a control law curve, Pe ¼ f ðNÞ; part of which exhibits an infinite slope (infinite derivative), as are the cases represented in Fig. 6 (without relief valve) and Fig. 8 (with relief valve), respectively. The maximum allowed slope depends on what the grid accepts in terms of power oscillation and on the inertia of the rotating parts, as will be shown in what follows. The dynamics of the rotor can be written as Pn ðtÞ 2 Pe ðtÞ ¼ IN

dN ; dt

Pn ¼ NLn ; Pe ¼ NLe ;

ð15Þ

where I is the inertia of the rotating parts (turbine and generator), Ln is the aerodynamic torque on the turbine rotor and Le is the electromagnetic torque on the generator rotor (bearing friction torque is ignored here). In the Pico plant it is I ¼ 595 kg m2 : From the grid viewpoint, the most unfavourable situation is expected to occur when the turbine torque sharply drops due to rotor blade stalling. Let us assume, in the worst scenario, that Ln has dropped to zero, and write 2Pe ¼ IN dN=dt; from which we obtain 2dPe = dt ¼ ðINÞ21 f ðNÞf 0 ðNÞ: We may specify 2dPe =dt # A; and consider the limiting case when 2dPe =dt ¼ A: We obtain an ordinary differential equation f ðNÞf 0 ðNÞ ¼ AIN: Let Nmax be the maximum value the rotational speed is allowed to take (for mechanical and/or electrical reasons). In order to avoid overspeeding, the control law should prescribe f ðNÞ ¼ Pmax ; N $ Nmax ; where Pmax is a value to be prescribed close to the peak power of the turbine at its maximum speed. The solution of the differential equation subject to boundary

80

A.F. de O. Falca˜o / Applied Ocean Research 24 (2002) 73–82

Fig. 10. Regulation curves for several values of AI £ 1023 (kW s21 kg m2), where A ¼ ldPe =dtl is the maximum value of the time-derivative of the electrical power allowed by the grid. The thick line is the optimal control algorithm for the Pico plant, without relief valve, for AI ¼ 50 £ 103 kW s21 kg m2 :

condition f ðNmax Þ ¼ Pmax is 2 2 N 2 Þ1=2 : Pe ¼ f ðNÞ ¼ ½P2max 2 AIðNmax

ð16Þ

This is the equation of a hyperbola in the ðN; Pe Þ plane, whose axes coincide with the axes of coordinates, passing through the points ðNmax ; Pmax Þ and ðN0 ; 0Þ; where N0 ¼ 2 ½Nmax 2 P2max ðAIÞ21 1=2 : Four regulation curves Pe ¼ f ðNÞ are shown in Fig. 10 for Pmax ¼ 500 kW; Nmax ¼ 157:1 rad s21 ; and AI £ 1023 ¼ 20; 30, 50, 100 kW s21 kg m2. As should be expected, the slope of the curve increases with the product AI. Obviously a larger slope allows more electrical energy to be produced; this may be achieved by increasing the rotational inertia I and/or the value of the electrical power time-derivative ldPe =dtl ¼ A allowed by the grid. We consider again the Pico plant, for which it is I ¼ 595 kg m2, and adopt the set of values Nmax ¼ 157:1 rad s21 ; Pmax ¼ 500 kW; and AI £ 1023 ¼ 50 kW s21 kgm2 (corresponding to A ¼ 84.0 kW s21). The control algorithm to be implemented in the plant’s PLC should be a combination of the curve Pe ¼ 1:224 £ 1025 N 3:1563 with the appropriate arc of hyperbola (thick line in Fig. 10).

4. Conclusions The stochastic modelling was found to be a fruitful approach to devise an optimal algorithm for the rotational speed control of an OWC plant equipped with a Wells turbine and to evaluate the average power output of the plant. The optimal control law may be expressed in terms of a simple relationship between the electromagnetic torque (to be applied on the generator rotor) and the rotational speed, and can be easily implemented into the generator power electronics and the plant’s programmable logical controller. This control is stable as was shown in Ref. [3] for the (slightly different) cube-law.

The electrical grid (especially in the case of a small isolated grid) may impose constraints on what concerns maximum allowable values for the time-derivative of the electrical power delivered to the grid. A simple algorithm, accounting for such restriction, was derived in order to complement the optimal control law. It is remarkable that the optimal rotational speed control algorithm is practically insensitive to the wave climate (as can be shown by the fact that all the points in Fig. 5, or in Fig. 7, are almost exactly aligned along a single curve, in spite of the fact that the corresponding points in Fig. 1 are widely scattered). For this reason it is to be presumed that the control algorithm would only be very weakly affected if the numerically computed excitation-volume-flow coefficient G(v ) were replaced by a more realistic curve accounting for losses due to real fluid effects in water (from the point of view of plant performance computation, such replacement would be equivalent to substituting the local wave climate by a less energetic one). On the other hand, the control algorithm is expected to depend strongly on the aerodynamic performance curves of the Wells turbine. So it is important to use realistic turbine performance curves (preferably from measurements in the prototype). We recall that several types of losses have been ignored, namely (1) real fluid, viscous losses in the water, and (2) all electrical losses (we have computed turbine, not electrical, power output). Besides, linearizing assumptions have been introduced, which are expected to lead to an overestimation of the produced energy: (3) linear, or small amplitude, water wave theory; (4) linearized isentropic relationship between air density and pressure in the chamber. The turbine aerodynamic performance was based on steady flow tests, and ignored hysteretic effects that could play a significant role in reciprocating flow (this is an area that requires further investigation).

Acknowledgements This work was partly funded by the JOULE Program of the European Commission under contracts No. JOR3-CT980282, JOR3-CT98-0312, and by the Portuguese Foundation for Science and Technology through IDMEC, Lisbon. The author would like to acknowledge Dr Ana Brito-Melo for supplying the results from numerical modelling on which the curves in Fig. 2 are based, and Dr Paulo Justino for helpful discussions.

Appendix A. Aerodynamic axial load on a Wells turbine rotor The air flow through the Wells turbine produces a reciprocating axial load, Fa, on the bearings, whose sign and

A.F. de O. Falca˜o / Applied Ocean Research 24 (2002) 73–82

magnitude depend on the instantaneous flow conditions. This force is difficult to measure either in model or at full scale, but can be evaluated approximately from theoretical considerations. Let lpl ¼ p0A 2 p0B ; where 0A and 0B denote conditions in the chamber and in the atmosphere, respectively, if p0A . p0B ; and vice versa if p0A , p0B : Assuming the flow to be incompressible ðr ø ra Þ; the axial Ð force on the rotor blades is ðp1 2 p2 ÞdA; where A ¼ pðD2 2 D2i Þ=4; D is the outer diameter, Di the inner diameter, and p1 and p2 are the pressure upstream, and the pressure downstream, of the rotor, respectively. On what follows we assume one-dimensional flow and write more simply ðp1 2 p2 ÞA: We consider first the case of a rotor without guide vanes. Neglecting the stagnation pressure losses between 0A and 1, we may write p1 ¼ p0A 2 ra Vx2 =2; where Vx ¼ mð _ ra AÞ21 is the axial component of the flow velocity. We assume that the whole kinetic energy is lost at the duct exit and write p2 ¼ p0B : It follows that p1 2 p2 ¼ lpl 2 ra Vx2 =2: We consider that the whole turbine rotor disc (including the rotor hub), r # D=2; is subject to the same pressure difference p1 2 p2 ; and so we find, in dimensionless form, the following expression for the axial load Fa lFa l p 1 ¼ lCl 2 F2 ; 4 2ð1 2 n2 Þ2 ra N 2 D4

ðA1Þ

where n ¼ Di =D; and C and F are the dimensionless coefficients of head and flow defined by Eqs. (3a) and (3b). We consider now one row of guide vanes on each side of the rotor, and assume that the aerodynamic design of the guide vane system follows the method developed in Ref. [13], based on two-dimensional incompressible irrotational flow through a triple cascade of blades (which approximately represents the non-stalled real flow). If the rotor blade profile is thin, then we may write a1 ¼ p 2 a2 ð0 , a1 , p=2Þ; where a1 ¼ arctanðVx =V1y Þ and a2 ¼ arctanðVx =V2y Þ are the flow angles at inlet to, and outlet from, the rotor, and V1y and V2y ð¼ 2V1y Þ are the corresponding values of the circumferential (swirl) component of the flow velocity. The difference Da ¼ a2 2 a1 ¼ p 2 2a1 is the angular deflection of the flow produced by the rotor blades. The theory shows that (as long as the flow remains irrotational, which obviously excludes stalling) the angles a1 and a2 are independent of the flow-rate m _ and rotational speed N (they depend only on the triple bladecascade geometry). Neglecting the aerodynamic losses at the guide vanes (where the flow velocity is much smaller that the blade speed), and keeping the same assumptions as for the single rotor geometry, we may write p1 ¼ p0A 2 ra

Vx2 csc2 a1 ; 2

p2 ¼ p0B 2 ra

Vx2 2 cot a1 ; 2

from which we find that Fa is still given by Eq. (A1). If severe stalling occurs at the rotor blades ðC . Ccr Þ (we denote the exit flow angle by a02 ), then the torque becomes small, and the same happens to the flow deflection

81

Da ¼ a02 2 a1 (we note that it is no longer a02 ¼ p 2 a1 ). In this case we assume approximately a02 ø a1 and V 02y ø V1y ¼ Vx tan a1 : The flow direction at the entrance to the second row of guide vanes is now far from its design value p 2 a1 and a substantial loss is expected to occur (due to large incidence), which we estimate to be equal to ra ðV2y 2 V 02y Þ2 =2 ¼ 2ra Vx2 cot2 a1 [14]. Combining these results, we obtain lFa l p 1 ¼ lCl 2 ½1 þ 4HðlCl 2 Ccr Þ 4 2ð1 2 n2 Þ2 ra N 2 D4 cot2 a1 F2 :

ðA2Þ

Here H(x ) is Heaviside’s unit step function, HðxÞ ¼ 0 for x , 0; HðxÞ ¼ 1 for x . 0: We recall that Eq. (A1) applies to the isolated rotor and Eq. (A2) to a turbine with guide vanes. We note that, if guide vanes are present, the angle a1 (directly related to guide vane geometry) is known to vary (in general) with the radial coordinate. The value of a1 that appears in Eq. (A2) should be understood as a radially averaged value. The above expressions, based on one-dimensional flow, should not be expected to provide very accurate evaluations of the axial load Fa. In the present paper it is assumed that the Wells turbine is a linear device, and so F ¼ K C; K depending on turbine shape, but not on its size, rotational speed or air density. We also assume that the bearing friction torque Ta is proportional to the axial load, i.e. Ta ¼ llFa l; l being a constant characterizing the bearing system (see Section 3.4). If, in addition, the (dimensionless) pressure-head C is a random variable with a Gaussian probability density function and standard deviation sC, then we may calculate the average value of Ta as a function of sC and find T
a ¼ lra N 2 D4

rffiffiffiffi p K 2 s2C sC 2 ð1 þ QÞ; 8 2ð1 2 n2 Þ2

ðA3Þ

where Q ¼ 0 for an isolated rotor, and   2 2x Q ¼ 4 cot2 a1 pffiffi e2x þ erfcx ; p

ðA4aÞ

C x ¼ pffiffi cr ; 2 sC

ðA4bÞ

for a rotor between guide vanes. If a relief valve is used, then, whenever lCl . Ccr ; lFa l is kept equal to its critical value lFa lcr (as given by Eq. (A2) for

A.F. de O. Falca˜o / Applied Ocean Research 24 (2002) 73–82

82

C ¼ Ccr ). We find rffiffiffiffi pffiffi 2 p T
a ¼ sC ð1 2 e2x þ x perfcxÞ 2 4 8 lra N D K 2 s2C 2ð1 2 n2 Þ2 "

2

2x e2x 1 þ ð2x 2 1Þerfcx 2 pffiffi p 2

2

# :

ðA5Þ

In Eqs. (A4a) and (A5), erfc is the complementary error function [15].

References [1] Falca˜o AF de O. The shoreline OWC wave power plant at the Azores. Proceedings of the Fourth European Wave Energy Conference, Aalborg, Denmark, 2000. p. 42–8. [2] Heath T, Whittaker TJT, Boake CB. The design, construction and operation of the LIMPET wave energy converter (Islay, Scotland). Proceedings of the Fourth European Wave Energy Conference, Aalborg, Denmark, 2000. p. 49–55. [3] Justino PAP, Falca˜o AF de O. Rotational speed control of an OWC wave power plant. Trans ASME: J Offshore Mech Arctic Engng 1999; 121:65–70.

[4] Falca˜o AF de O, Rodrigues RJA. Stochastic modelling of OWC wave power plant performance. Appl Ocean Res, this issue (PII: S01411187(02)00022-6). [5] Kinsman B. Wind waves. Englewood Cliffs, NJ: Prentice-Hall; 1965. [6] Sarpkaya T, Isaacson M. Mechanics of wave forces on offshore structures. New York: Van Nostrand Reinhold; 1981. [7] Goda Y. Random seas and design of maritime structures, 2nd ed. Singapore: World Scientific; 2000. [8] Dixon SL. Fluid mechanics and thermodynamics of turbomachinery, 4th ed. Boston: Butterworth –Heinmann; 1998. [9] Falnes J. Ocean waves and oscillating systems. Cambridge: Cambridge University Press; 2002. [10] Leon-Garcia A. Probability and random processes for electrical engineering, 2nd ed. Reading: Addison-Wesley; 1994. [11] Brito-Melo A, Hofmann T, Sarmento AJNA, Cle´ment AH, Delhommeau G. Numerical modelling of OWC-shoreline devices including the effect of surrounding coastline and non-flat bottom. Int J Offshore Polar Engng 2001;11:147–54. [12] Falca˜o AF de O, Justino PAP. OWC wave energy devices with air flow control. Ocean Engng 1999;26:1275–95. [13] Gato LMC, Falca˜o AF de O. Performance of Wells turbine with double row of guide vanes. Int J Jpn Soc Mech Engrs, Series II 1990; 33:265–71. [14] Pfleiderer C, Petermann H. Stro¨mungsmaschinen, 6th ed. Berlin: Springer; 1991. [15] Abramowitz M, Stegun IA. Handbook of mathematical functions. New York: Dover; 1965.