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 coefﬁcients 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 be written as P ¼ W w 2 L ; where Ww is the power absorbed from the waves and L are the losses that occur in the energy conversion chain. The most signiﬁcant losses are: (i) L1 due to real (viscous) ﬂuid 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 variable- E-mail address: falcao@hidro1.ist.utl.pt (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 efﬁciency depends on ﬂow-rate and rotational speed. On the other hand, variations in the turbine rotational speed affect the ﬂow-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 ﬂow is assumed incompressible and the effects due to Reynolds number and Mach number variation neglected, the instantaneous value of the aerodynamic efﬁciency of a given turbine is known to be a function only of Tt=N2; 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 £ N2 to be implemented in the programmable logic controller (PLC) of the plant. This approach aimed at keeping the timeaveraged value of the turbine aerodynamic efﬁciency 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 efﬁciency 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 signiﬁcant, but are difﬁcult 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 ﬂow 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 brieﬂy 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 coefﬁcients 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 ﬂow-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 ﬁxed 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 ! f ðzÞ ¼ pﬃﬃ1ﬃﬃ exp 2psz z2 2 2s2z ; ð1Þ where s2z is the variance, and sz is the standard deviation, of z. We assume here one-directional incident waves (multi- directional 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) deﬁned in the interval ð21; 1Þ and is an even function. It is related to the more frequently used one-sided spectral power density Sðz1ÞðvÞ by Sðz1ÞðvÞ ¼ 2SzðvÞ if v . 0; Sðz1Þð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 D2 ; ð3aÞ F¼ m_ raND3 ; ð3bÞ P¼ Pt ra N 3 D5 : ð3cÞ Here C is non-dimensional pressure-head, F is nondimensional ﬂow-rate, P is non-dimensional turbine power output, ra is atmospheric air density, m_ is mass ﬂow-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 75 by F ¼ KC; ð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Þ 0 Here G(v ) is the excitation-volume-ﬂow coefﬁcient (as deﬁned in Ref. [9]), and L ¼ KD þ B þ i vV0 21 þC ; ð7Þ raN gpa where pa is the atmospheric pressure, g ¼ cp=cv is the speciﬁc heat ratio for air, V0 is the undisturbed value of the air chamber volume, and the hydrodynamic coefﬁcients 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 P t ¼ rpaﬃNﬃﬃ3ﬃD5 2psp ð1 exp 21 2 p2 2s2p ! fP p ra N 2 D2 dp; ð8Þ or in dimensionless form P ¼ pﬃﬃ1ﬃﬃ 2psC ð1 exp 21 ! C2 2 2s2C fPðCÞdC; ð9Þ where P ¼ P t ra N 3 D5 ; sC ¼ sp ra N 2 D2 : ð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 ﬁnd the optimal relationship, for each sea state, between the electromagnetic torque Te (or the corresponding power Pe ¼ TeN) 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 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 deﬁnitions (10a) and (10b), we easily obtain ! 1 dP t ¼ 3P þ D3sp dN sC N dsp 2 2 sp dN dP : dsC ð11Þ The condition for maximum turbine power output implies dP t=dN ¼ 0; and so, from Eq. (11), it follows that sC P dP dsC ¼ 22 3 N dsp : ð12Þ sp dN 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 Pp: From Eq. (10a) we ﬁnd P t ¼ raD5PpN3: ð13Þ We note that Pp 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 £ N3 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 ﬂow 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 signiﬁcant 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 coefﬁcients for the Pico plant: radiation conductance B (solid line), radiation susceptance C (dashed line) and excitation-volumeﬂow coefﬁcient 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 ﬂow (chain line). whose area is proportional to the sea state contribution to the total annual energy ﬂux 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:5Hs2Te24v25 expð21054Te24v24Þ; ð14Þ where Hs is the signiﬁcant wave height and Te is the energy period. Curves for the hydrodynamic coefﬁcients 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 reﬂects 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 ¼ KC; 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 inﬁnity, 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 77 Fig. 4. Time-averaged bearing friction power loss NTa; 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 ﬂow 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 ﬂow 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 ¼ KC; is still assumed between the dimensionless excess pressure in the chamber C and the dimensionless volume ﬂow rate F through the turbine.) 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-volumeﬂow coefﬁcient G depends only on the chamber geometry, whereas L depends also on KX (where X ¼ DðraNÞ21) and 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ðraNÞ21Þ: Since the dimensionless timeaveraged power output of the turbine P is a known function of sC (Fig. 3), it can be regarded (for ﬁxed ra and D ) as a function of N. The time-averaged turbine power can be calculated from P t ¼ raN3D5P: 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 ﬂow 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 ﬂow-rate through the turbine; this is difﬁcult to measure in practice, but can be approximately computed from theoretical ﬂow 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=lFal ¼ 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 NTa (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 NTa as shown in Fig. 4 (compare the solid and broken lines). 3.5. Numerical results 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 ﬂow 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 £ 1025N3:1616; with P n in kW and N in rad s21. 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 ﬂow 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 £ 1025N3: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 signiﬁcantly. It is to be expected that the points representing the optimal plant performance under this larger set of Fig. 7. As in Fig. 5, for a plant with ﬂow control by relief valve. The curve’s equation is Pe ¼ P n ¼ 1:093 £ 1025N3:2361: 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 speciﬁed 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 modiﬁed, 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 ﬁnd, 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 £ 1025N3:1563 ( ¼ 1.093 £ 1025N 3.2361). Dashed lines represent simpliﬁed laws Ppe ¼ 3:344 £ 1025N3 ( ¼ 5.054 £ 1025N 3). the case when a controlled relief valve prevents the air ﬂow rate through the turbine from exceeding the stall-free limit. The new control equation is now Pe ¼ P n ¼ 1:093 £ 1025N3: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 £ 1025N3:1563; shown in Fig. 5, with the results from the simpliﬁed control law given by Eq. (13). We ﬁnd spC ¼ 0:03253 and Pp ¼ 0:0004157; from which we obtain PpeðNÞ ¼ 3:344 £ 1025N3: If a relief valve is operating, the corresponding control laws are Pe ¼ P n ¼ 1:093 £ 1025N3:2361 and PpeðNÞ ¼ 5:054 £ 1025N3: 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 ﬂow is controlled by a relief valve. The dashed (thick and thin) lines represent the simpliﬁed control laws. It may be seen that the N 3-laws (dashed lines as compared with solid lines) specify signiﬁcantly 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 X44 P n;annual ¼ fjP 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 ﬂow 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 ﬂywheel 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 inﬁnite slope (inﬁnite 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Þ21f ð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 Pe ¼ f ðNÞ ¼ ½P2max 2 AIðNm2 ax 2 N2Þ1=2: ð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 ¼ ½Nm2 ax 2 P2maxðAIÞ211=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 s21kgm2 (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 £ 1025N3: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-ﬂow coefﬁcient G(v ) were replaced by a more realistic curve accounting for losses due to real ﬂuid 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 ﬂuid, 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 ﬂow tests, and ignored hysteretic effects that could play a signiﬁcant role in reciprocating ﬂow (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 ﬂow 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 81 magnitude depend on the instantaneous ﬂow conditions. This force is difﬁcult 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 ﬂow to be incompÐressible ð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 ﬂow and write more simply ðp1 2 p2ÞA: We consider ﬁrst the case of a rotor without guide vanes. Neglecting the stagnation pressure losses between 0A and 1, we may write p1 ¼ p0A 2 raVx2=2; where Vx ¼ m_ ðraAÞ21 is the axial component of the ﬂow 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 raVx2=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 ﬁnd, in dimensionless form, the following expression for the axial load Fa lFal ra N 2 D4 ¼ p lCl 2 4 2ð1 1 2 n2Þ2 F2; ðA1Þ where n ¼ Di=D; and C and F are the dimensionless coefﬁcients of head and ﬂow deﬁned 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 ﬂow through a triple cascade of blades (which approxi- mately represents the non-stalled real ﬂow). If the rotor blade proﬁle 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 ﬂow 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 ﬂow velocity. The difference Da ¼ a2 2 a1 ¼ p 2 2a1 is the angular deﬂection of the ﬂow produced by the rotor blades. The theory shows that (as long as the ﬂow remains irrotational, which obviously excludes stalling) the angles a1 and a2 are independent of the ﬂow-rate m_ and rotational speed N (they depend only on the triple blade- cascade geometry). Neglecting the aerodynamic losses at the guide vanes (where the ﬂow 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 2 csc2a1; p2 ¼ p0B 2 ra Vx2 2 cot2a1; from which we ﬁnd that Fa is still given by Eq. (A1). If severe stalling occurs at the rotor blades ðC . CcrÞ (we denote the exit ﬂow angle by a02), then the torque becomes small, and the same happens to the ﬂow deﬂection 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 0 2y ø V1y ¼ Vx tan a1: The ﬂow 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 V02yÞ2=2 ¼ 2raVx2 cot2a1 [14]. Combining these results, we obtain lFal ra N 2 D4 ¼ p lCl 2 4 2ð1 1 2 n2Þ2 ½1 þ 4HðlCl 2 CcrÞ cot2a1F2: ð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 ﬂow, 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 ¼ KC; 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 ¼ llFal; 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 ﬁnd T a lra N 2 D4 rﬃﬃﬃﬃ ¼ p 8 sC 2 K 2 s2C 2ð1 2 n2Þ2 ð1 þ QÞ; ðA3Þ where Q ¼ 0 for an isolated rotor, and Q ¼ 4 cot2a1 p2ﬃxﬃ e2x2 þ erfcx ; p ðA4aÞ x ¼ pCﬃﬃ cr ; 2sC ðA4bÞ for a rotor between guide vanes. If a relief valve is used, then, whenever lCl . Ccr; lFal is kept equal to its critical value lFalcr (as given by Eq. (A2) for 82 A.F. de O. Falca˜o / Applied Ocean Research 24 (2002) 73–82 C ¼ Ccr). We ﬁnd T a lra N 2 D4 ¼ rﬃﬃﬃﬃ p 8 sC ð1 2 e2x2 þ pﬃﬃ x perfcxÞ 2 K 2 s2C 2ð1 2 n2Þ2 " # Â 1 þ ð2x2 2 1Þerfcx 2 2xpeﬃ2ﬃ x2 : p ð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 Scientiﬁc; 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: Cam- bridge 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, Delhom- meau G. Numerical modelling of OWC-shoreline devices including the effect of surrounding coastline and non-ﬂat bottom. Int J Offshore Polar Engng 2001;11:147–54. [12] Falca˜o AF de O, Justino PAP. OWC wave energy devices with air ﬂow 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] Pﬂeiderer C, Petermann H. Stro¨mungsmaschinen, 6th ed. Berlin: Springer; 1991. [15] Abramowitz M, Stegun IA. Handbook of mathematical functions. New York: Dover; 1965.

## 评论