SOME QUESTIONS OF WIND-ENERGY CALCULATIONS

V. M. Osipov

(Presented by Professor, Doctor of Technical Sciences
I. D. Kutyavin)

In recent years, in our country and abroad, serious successes have been achieved in the technology of using wind energy. New, more advanced designs of wind engines with high technical and economic indicators have been created. The studies carried out and the available operating experience quite clearly show the economic expediency of using wind energy for purposes of rural electrification. The broad use of wind energy to increase the power supply of socialist agriculture through the operation of wind power stations (WPS) in local systems, as well as in parallel with the powerful networks of state power systems, poses a number of new tasks for wind-power specialists. These tasks are conditioned by the special nature of wind energy, the use of which in an energy system generally creates special features with respect to its operating regime and affects its technical and economic indicators. The study of these questions is inevitable with the broad use of WPS for electrification purposes. The existing methodology of wind-energy calculations is not capable of solving the tasks posed, owing to the weakness of its theoretical basis. Moreover, according to the existing methodology, the quantities determined are obtained as the result of summing several terms, each of which corresponds to an integral value of wind speed, since the statistical set for wind recurrence is usually given in the form of a broken line, and the principal characteristic of the wind engine has no analytical expression. Calculations are usually carried out in the form of tables, and the variant principle is widely used. These circumstances make such calculations very labor-intensive, especially when they must be performed for several values of the mean annual wind speed or other passport data of the WPS; moreover, obviously, analysis of the influence of various factors is extremely difficult and can be carried out only by comparing different calculation variants. Thus, the existing methodology of wind-energy calculations is theoretically very imperfect, and practically inconvenient because of its laboriousness.

In light of the foregoing, it appears to us an urgent task to change the existing methodology with the aim of simplifying it as much as possible while preserving the necessary accuracy, and also with the aim of applying it to the solution of more complex questions of wind utilization.

In the present work an attempt is made in this direction. In doing so, we proceed from the following premises:

1. The wind speed, which determines its energy possibilities, is a continuous random variable capable of varying over a very wide range of values. Consequently, the power of the wind flow is also a continuous random variable; this, as is known, is precisely the special feature of wind energy.

1. Wind, as an energy carrier, by virtue of what has been said, is subject to the regularities of random variables and, consequently, probability theory— the branch of mathematics that studies these regularities—must serve as the theoretical foundation of wind-power engineering.

Basic Concepts

A random variable is considered to be defined if some function is known that has the property that the product \(f(v)\Delta v\) is equal to the probability that the random variable \(v\) will fall in the interval between \(v\) and \(v+\Delta v\) [1, 2]. On the other hand, according to the mathematical (statistical) concept of probability [1], the probability that the wind speed \(v\) will lie within the interval \(v+\Delta v\) is equal to the ratio \(\frac{\Delta \tau}{T}\), where \(\Delta \tau\) is the number of hours in a year during which all speeds of the indicated interval will be observed; \(T\) is the annual number of hours.

Thus, we have:

\[ \frac{\Delta \tau}{T}=f(v)\Delta v, \tag{1} \]

or, dividing both sides of the equality by \(\Delta v\) and passing to the limit, we obtain:

\[ \frac{1}{T}\cdot \frac{d\tau}{dv}=f(v). \tag{2} \]

The function \(f(v)\), defined in this way, is called the probability-density distribution function of the random variable \(v\).

In wind-power engineering, use is made of a certain value of probability density averaged over the interval \(\Delta v\), with \(\Delta v=1\) being chosen; then the density averaged over the interval

\[ \frac{\Delta \tau}{T\Delta v}=\frac{\Delta \tau}{T}, \]

i.e., turns out to be numerically equal to the probability that the quantity \(v\) will fall in the given interval \(\Delta v=1\); this latter is also called the frequency.

Thus, frequency is nothing other than the probability density averaged over the interval \(\Delta v=1\), numerically equal to the probability of the inequality \(B(v\leq \varepsilon \leq v+1)\), i.e., the probability that the quantity \(v\) will fall in the interval between \(v\) and \(v+1\).

Having the distribution function, it is easy to obtain the basic numerical characteristics of the random variable—the moments. The moment of order \(k\) is equal to the mathematical expectation of the \(k\)-th order [1, 2], i.e.,

\[ M_k(v)=\int_a^b v^k f(v)\,dv. \tag{3} \]

From the definition of the distribution function it follows that the moment of order zero is equal to the probability of the inequality

\[ M_0(v)=B(a\leq v\leq b)=\int_a^b f(v)\,dv. \]

If the integration covers all observed speeds, then, obviously:

\[ B(a\leq v\leq b)=\int_a^b f(v)\,dv=1, \]

i.e., the area under the distribution curve is equal to unity.

A random variable can be defined by using a function other than the probability-density distribution function \(f(v)\). The graph of this function is widely used in power engineering for estimating the load over an annual interval and has received the name of a duration curve. The latter gives a clear idea of how many hours in a year the values of the speed greater and less than a given value \(v_1\) will be observed (Fig. 1).

The quantity \(t_{v_1}\) denotes the annual duration of speeds greater than \(v_1\), and \(T-t_{v_1}\) the duration of speeds less than \(v_1\). In relative units (i.e., \(t_{v_1}/T\) and \((T-t_{v_1})/T\)) this will mean the probability of encountering a speed, respectively, greater than \(v_1\) and less than \(v_1\). Thus we may write the obvious equality:

Fig. 1.

\[ \frac{T-t_{v_1}}{T}=B(0\le v\le v_1). \]

On the other hand,

\[ B(0\le v\le v_1)=\int_0^{v_1} f(v)\,dv \]

and, consequently,

\[ \frac{T-t_{v_1}}{T}=\int_0^{v_1} f(v)\,dv. \]

In general we may write:

\[ \frac{T-t}{T}=\int_0^v f(v)\,dv . \tag{4} \]

Differentiating with respect to the variable upper limit, we obtain:

\[ \frac{d(T-t)}{T\,dv}=-\frac{dt}{T\,dv}=f(v). \tag{5} \]

The last equality establishes the connection between the equation of the duration curve and the probability-density distribution function. It proves very useful for our purposes.

From (4) we shall have:

\[ t=t(v)=T\left[1-\int_0^v f(v)\,dv\right]. \]

Let us denote:

\[ F(v)=\int_0^v f(v)\,dv, \tag{6} \]

then we obtain:

\[ t=T[1-F(v)]. \tag{7} \]

The function \(F(v)\) is called the integral distribution function \([1,2]\). It evidently determines the probability that our random variable will be less than a given value \(v\), i.e. is equal to the probability of the inequality \(B(\varepsilon<v)\).

As \(v\) increases, the probability of this inequality will increase and, consequently, \(F(v)\) is a nondecreasing function \([2]\).

It follows from (7) that the ratio

\[ \frac{t}{T}=\varphi(v)=1-F(v) \tag{8} \]

means the probability of the opposite event [i.e., the probability of the inequality \(B(\varepsilon>v)\)] and is, consequently, the integral distribution function of the probability of the opposite event.

The function \(\varphi(v)\), in contrast to \(F(v)\), will be a decreasing function. We shall call it the decreasing integral distribution function. Bearing in mind that \(dF(v)=f(v)\,dv\), we can write the moment of order \(k\) in the form of a Stieltjes integral:

\[ M_k=\int_a^b v^k\,dF(v). \]

Noting further that \(T\,dF(v)=-dt\) and that to the upper limit there corresponds some quantity \(t_1\), and to the lower limit the quantity \(t_2\), with \(b>a\), \(t_1<t_2\), we obtain:

\[ M_k=-\frac{1}{T}\int_{t_2}^{t_1} v^k\,dt =\frac{1}{T}\int_{t_1}^{t_2} v^k\,dt . \tag{9} \]

The last integral represents the area under a certain curve obtained as a result of raising the ordinates of the decreasing integral distribution curve \(\dfrac{t}{T}=t_*= \varphi(v)\) (i.e., the ordinates of the duration graph) to the \(k\)-th power. This is the geometrical meaning of the moment of order \(k\).

If the random variable \(v\), having probability density \(f(v)\), enters as the argument into some function \(P=\psi(v)\), then the latter, obviously, will also be a random variable with some density \(f_1(P)\). The relation between the densities is expressed by the equality:

\[ f(v)\,dv=f_1(P)\,dP . \tag{10} \]

It establishes the fact that the corresponding area elements under the distribution curves will be equal [3]. From (10) we shall have:

\[ f_1(P)= \left. \frac{f(v)}{\dfrac{dP}{dv}} \right|_{v=\psi^{-1}(P)} . \tag{11} \]

On the same basis one may write:

\[ \int_0^v f(v)\,dv = \int_{P_0=\psi(0)}^{P=\psi(v)} f_1(P)\,dP \tag{12} \]

and, consequently,

\[ t=t(v)=T\left[1-\int_0^v f(v)\,dv\right] = T\left[1-\int_{P_0}^{P} f_1(P)\,dP\right]=t(P), \]

i.e.

\[ t(v)=t(P). \]

Thus, the duration graphs of two functionally related random variables differ only in their ordinates; the abscissas of the corresponding points are identical. To each value \(v\) there will correspond perfectly definite values of \(P\) and \(t\); consequently,

\[ \begin{gathered} P=\psi(v);\\ t=T\left[1-\int_0^v f(v)\,dv\right]. \end{gathered} \tag{13} \]

is the equation of the duration curve in parametric form for the random variable \(P\) (the equation of the decreasing integral distribution curve). The role of the parameter is played by the quantity \(v\). Taking into account that

\[ \int_{0}^{v_m} f(v)\,dv = \int_{0}^{v} f(v)\,dv + \int_{v}^{v_m} f(v)\,dv = 1, \]

where \(v_m\) is the maximum value of the wind speed, the dependence \(t(v)\) can be represented in the form:

\[ t = T \int_{v}^{v_m} f(v)\,dv . \tag{14} \]

Fig. 2.

In accordance with this, one may assert that

\[ \left. \begin{aligned} P_v &= c_0 D^2 v^3 \xi \ [\mathrm{kW}];\\ t &= T \int_{v}^{v_m} f(v)\,dv \end{aligned} \right\} \tag{15} \]

is the equation of the duration curve for the power of a wind motor in parametric form. The following notations have been adopted here:

\(c_0\) — a constant determined by the density of the air and the selected system of units (usually \(c_0 = 0.481 \cdot 10^{-3}\) is taken);

\(D\) — diameter of the wind wheel, m;

\(v\) — wind speed, m/sec;

\(\xi\) — coefficient of utilization of wind energy (in the general case, a function of the circumferential speed of rotation of the wind wheel and of the wind speed).

The work performed by the wind motor during the year is, evidently, equal to the area under the curve \(P_v(t)\) (Fig. 2), bounded between zero and the design power, i.e.

\[ A = \int_{0}^{P_{vp}} t\,dP_v \ [\mathrm{kW \cdot h}]. \tag{16} \]

Integrating by parts and interchanging the limits of integration, we obtain:

\[ A = P_{vp} t_p + \int_{t_p}^{t_0} P_v\,dt \ [\mathrm{kW \cdot h}], \tag{17} \]

where \(P_{vp}\) is the design power of the wind motor;

\(t_p\) — number of hours of operation of the installation at design power;

\(t_0\) — total number of hours of operation.

The annual output may be expressed through functions of the probability-density distribution. Having the equation of the duration curve for the power of the wind motor, it is easy to obtain the moment of any order. The incomplete moment of the \(k\)-th order, taken within the limits from zero to \(P_{vp}\), as follows from the geometrical meaning of the moment, may be determined by the formula

\[ M_k = P_{vp}^{k}\frac{t_p}{T} + \frac{1}{T}\int_{t_p}^{t_0} P_v^{k}\,dt . \]

Integrating by parts, we easily obtain:

\[ M_k=\frac{1}{T}\int_{0}^{P_{vr}} t\,d\left(P_v^k\right). \tag{18} \]

These are some of the general theoretical relations that must be laid as the basis for wind-power calculations.

Distribution functions of wind speed

For carrying out specific calculations it is necessary to know the distribution function of the probability density (frequency) of wind speed. Such a function can be determined by processing statistical data from observations of wind speed at some station or group of stations. In Soviet wind power engineering, two different distribution functions, obtained by processing the statistical material of a number of meteorological stations, have mainly become widespread. One of these functions was proposed by M. M. Pomortsev [6, 7] and is a normal law; the other—by Gullen—in the form:

\[ f(v)_z=k_1e^{-3av}\sin(av). \tag{19} \]

For some regions of our country one law is confirmed, for others—the other [6, 7]. The recommendations on this matter available in the literature are apparently insufficiently substantiated. According to M. M. Pomortsev, the density obeys the normal law, which plays an enormous role in probability theory:

\[ f(v)_n=k_2 e^{-\frac{(v-b)^2}{2\sigma^2}}. \tag{20} \]

The coefficients \(k_1\) and \(k_2\) are found from the condition that the area under the curve, encompassing all possible values of the random variable, be equal to unity. In the case of the normal law, the area must be calculated within the limits from \(-\infty\) to \(+\infty\); then

\[ k_2=\frac{1}{\sigma\sqrt{2\pi}}. \tag{21} \]

The quantity \(b\) in formula (20), called the center of the distribution, is equal to the most probable value of the wind speed. If \(k_2\) is chosen in accordance with (21), i.e., if the negative region is taken into account, then \(b=v_{sr}\). In reality, wind speed cannot take negative values; therefore \(b\) differs somewhat from \(v_{sr}\), but this difference is small owing to the insignificant influence of the negative region. This influence decreases rapidly with increasing \(v_{sr}\) [7]. The quantity \(\sigma\) is called the root-mean-square deviation and is calculated by the formula

\[ \sigma=\sqrt{v_{sr.kv}^{\,2}-v_{sr}^{\,2}}. \tag{22} \]

The square of the root-mean-square deviation is called the variance, i.e.,

\[ D(v)=\sigma^2. \tag{23} \]

Thus, the normal law is determined by two parameters: \(v_{sr}\) and \(\sigma\). However, in our case, as calculations show, the root-mean-square deviation \(\sigma\) is a function of the annual mean speed \(v_{sr}\). A graph of this dependence is given in the drawing (Fig. 3).

Above we assumed that wind speed can take infinitely large values. In reality, the maximum possible value of wind speed is a finite quantity; therefore, assuming,

that \(f(\infty)=0\), whereas the density becomes zero at a finite value of \(v\), we thereby allow a certain error. It is not difficult, however, to verify that this error in all practical cases will be a very small quantity.

The Gullen density becomes zero at a wind speed equal to zero, and at \(av=\pi\) (in general at \(av=n\pi\), where \(n=0,1,2,\ldots\)). For \(av>\pi\), \(f(v)_2<0\), which makes no sense, since a probability density is an essentially positive quantity. This occurs because the Gullen curve, like the Pearson curves [3], has no other justification than that these curves well reflect one or another statistical regularity [2, 3]. If we put \(av_m=\pi\), then

Fig. 3.

\[ a=\frac{\pi}{v_m}. \]

For \(av>\pi\)

\[ |f(v)_2|<k_1\cdot 0.000083. \]

It is clear that such a quantity may be neglected, i.e., one may assume

\[ f(v)_2\big|_{av>\pi}=0, \]

which is equivalent to admitting \(v_m=\infty\).

Now the parameter \(a\) can be found from the condition

\[ k_1\int_0^\infty v\exp(-3av)\sin(av)\,dv=v_{cp}, \]

and the coefficient \(k_1\) from the condition

\[ k_1\int_0^\infty \exp(-3av)\sin(av)\,dv=1. \]

This gives

\[ a^2=\frac{3k_1}{50v_{cp}} \]

and

\[ k_1=10a. \]

Solving jointly, we obtain:

\[ a=\frac{3}{5v_{cp}} \tag{24} \]

and

\[ k_1=10a=\frac{6}{v_{cp}}. \tag{25} \]

The obtained values of \(a\) and \(k_1\) practically coincide with Emtsev’s data, found by him in a somewhat different way [6].

Let us now find the equations of the cumulative curves for our two distributions. For the Gullen distribution this presents no difficulty,

\[ F(v)_2=k_1\int_0^v \exp(-3av)\sin(av)\,dv=\{1-\exp(-3av)[3\sin(av)+\cos(av)]\}. \]

Substituting into (7), we obtain:

\[ t=T\exp(-3a\upsilon)[3\sin(a\upsilon)+\cos(a\upsilon)]. \tag{26} \]

For various annual-average velocities we obtain a family of curves

\[ t_*=\frac{t}{T}=t_*(\upsilon); \]

it is more convenient, however, to take the quantity \(a\upsilon\) as the argument; then, instead of a family, we obtain a single curve, shown in Fig. 4.

The integral distribution \(F(\upsilon)\) in the case of the normal law has the form:

\[ F(\upsilon)_n=\frac{1}{\sigma\sqrt{2\pi}} \int_{-\infty}^{\upsilon} \exp\left[-\frac{(\upsilon-\upsilon_{cp})^2}{2\sigma^2}\right]\,d\upsilon. \tag{27} \]

Fig. 4.

This integral, because of its transcendental nature, cannot be taken in elementary functions. By the substitution \(y=\dfrac{\upsilon-\upsilon_{cp}}{\sigma}\) we can reduce it to the form:

\[ F\left(\frac{\upsilon-\upsilon_{cp}}{\sigma}\right) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{(\upsilon-\upsilon_{cp}):\sigma} e^{-\frac{y^2}{2}}\,dy = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{0} e^{-\frac{y^2}{2}}\,dy + \frac{1}{\sqrt{2\pi}} \int_{0}^{(\upsilon-\upsilon_{cp}):\sigma} e^{-\frac{y^2}{2}}\,dy. \]

Further, bearing in mind that

\[ \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{0} e^{-\frac{y^2}{2}}\,dy = \frac{1}{2}, \]

we obtain:

\[ F\left(\frac{\upsilon-\upsilon_{cp}}{\sigma}\right) = \frac{1}{2} \left[ 1+\frac{2}{\sqrt{2\pi}} \int_{0}^{(\upsilon-\upsilon_{cp}):\sigma} e^{-\frac{y^2}{2}}\,dy \right]. \tag{28} \]

This equality is valid when \(\dfrac{\upsilon-\upsilon_{cp}}{\sigma}>0\) \((\upsilon>\upsilon_{cp})\). In the case when \(\upsilon<\upsilon_{cp}\), we shall have [2]:

\[ F(-y)=1-F(y). \]

The integral

\[ \frac{2}{\sqrt{2\pi}} \int_{0}^{(\upsilon-\upsilon_{cp}):\sigma} e^{-\frac{y^2}{2}}\,dy = \Phi\left(\frac{\upsilon-\upsilon_{cp}}{\sigma}\right) \]

is called the probability integral; by the substitution \(s=\dfrac{y}{\sqrt{2}}\) it is reduced to the Laplace–Kramp function [1, 4, 5]. For the probability integral there exist detailed tables, with the aid of which \(F\left(\dfrac{\upsilon-\upsilon_{cp}}{\sigma}\right)\) can

can be constructed. The equation of the duration curves takes the form:

\[ \begin{gathered} t_*=\frac{1}{2}\left[1+\Phi\left(\frac{v-v_{\mathrm{av}}}{\sigma}\right)\right]\quad \text{for } v>v_{\mathrm{av}},\\ t_*=\frac{1}{2}\left[1-\Phi\left(\frac{v-v_{\mathrm{av}}}{\sigma}\right)\right]\quad \text{for } 0<v<v_{\mathrm{av}}. \end{gathered} \tag{29} \]

For \(v=v_{\mathrm{av}}\), \(\Phi(0)=0\) and \(t_*=\dfrac{1}{2}\).

In Fig. 5, according to equations (29), a family of duration curves has been constructed for various mean annual wind speeds (the lower curve for \(v_{\mathrm{av}}=5\ \mathrm{m/sec}\), the upper one for \(v_{\mathrm{av}}=10\ \mathrm{m/sec}\)).

Fig. 5.

Comparing the duration curves for the theoretical distributions with one another, one can see that the Gullen curve is characterized by the absence of a calm period and by an increase in duration in the region of high speeds in comparison with the normal distribution. Having these curves, it is easy to obtain duration graphs for the power of a wind turbine or a wind power station. The latter exhaustively characterize the annual output of a wind power station and the conditions of its operation. The relations and formulas obtained make it possible to solve various questions of the utilization of wind energy. Unfortunately, for lack of space the author is deprived of the opportunity to illustrate the application of the formulas and dependences obtained to the solution of specific problems of wind power engineering, above all such as the determination of the annual output of a wind power station in the most advantageous aerodynamic regime, and of other quantities characterizing the operating conditions of a wind power station.

Conclusions

1. The existing methodology of wind-power calculations is theoretically very imperfect and practically inconvenient because of its laboriousness.

2. Wind as an energy carrier obeys the laws of random variables and, consequently, probability theory, which studies these laws, must serve as the theoretical basis for wind-power calculations.

3. Wind recurrence is nothing other than the average, over the interval \(\Delta v=1\), probability density, numerically equal to the probability of the inequality \(B(v<\varepsilon<v+1)\).

4. It is expedient to proceed from the decreasing integral distribution function used in power engineering under the name of the duration curve. This creates clarity in the use of the methods of probability theory.

5. The distribution function of the power of a wind turbine as a random variable can easily be expressed through the distribution function of wind speed; in this case it is especially simple to obtain the equation of the duration curve.

6. A general expression for determining the annual output can be obtained on the basis of the duration curve.

1. The duration curve of the power of a wind engine or wind power station (WPS) exhaustively characterizes the efficiency and operating conditions of a wind-power installation; therefore, it should be taken as the basis when solving various questions concerning the use of wind energy.

References

1. Gnedenko B. V., Course in Probability Theory, Gostekhizdat, Moscow, 1954.

2. Boes G. P., Probability Theory, Gostekhizdat, Moscow—Leningrad, 1950.

3. Thornton Fry, Probability Theory for Engineers, translated from the English, edited by A. Ya. Khinchin, ONTI, Moscow—Leningrad, 1934.

4. Posse K., Course in Integral Calculus, GONTI, 1938.

5. Ryzhik I. M., Tables of Integrals, Sums, Series, and Products, Gostekhizdat, Moscow—Leningrad, 1952.

6. Fateev E. M., Wind Engines and Wind Installations, OGIZ, Selkhozgiz, Moscow, 1948.

7. Rosentul S., Design Parameters and Technical Operating Characteristics of Wind Installations, dissertation, Moscow, 1940.