350 lines
44 KiB
Markdown
Raw Normal View History

2025-09-08 17:02:40 +08:00
# SYMBOLIC COMPUTING AS A TOOL IN WIND TURBINE DYNAMICS
A. D. GARRADt AND D. C. QUARTON
Wind Energy Group, Taylor Woodrow Construction, 345 Ruislip Road, Southall UB1 2QX, England (Received 26 October 1984, and in revised form 25 September 1985)
The problem of deriving equations of motion that describe the coupled tower-rotor system of a wind turbine is discussed. The mathematical formulation of these equations is described first in a form suitable for manual derivation and then as a step by step process suitable for automation. Reasons for and experience in using a symbolic computing system to undertake this work for large models are described. The general approach is illustrated by means of a very simple model. Finally, some results of the stability analysis of the simple example and an eleven degree of freedom three-bladed model are presented together with a set of measured and predicted blade loads for a $250\,\mathbf{kW}$ wind turbine.
# 1. INTRODUCTION
We recently encountered a problem of algebraic manipulation in wind turbine dynamics which was conceptually very straightforward but immensely tedious to perform. We required a solution to this problem and were interested in finding a means of performing it automatically, first as a check on the manual calculation, and then as a method in its own right. The power and possible application'of computers as “"equation crunchers" rather'than "number crunchers" was a revelation to us. In this paper we shall not seek to break any new ground, but rather describe the analysis which we conducted as an illustration of the possible applications of symbolic or algebraic computing as a tool in engineering analysis. Our particular problem was concerned with formulating the equations of motion of a wind turbine system. This problem will be used here as an example, but it will become obvious that many problems in engineering dynamics are of a similar nature.
Symbolic computing systems have been in existence for some 20 years. References were made to this application of digital computers as early as 1953 [1, 2]. The systems seem to have been developed for high energy physics applications and, indeed, the system used by us (REDUCE) bears definite signs of its origins in this respect. Their application to engineering in general and engineering mechanics in particular seems to have been very limited. Some work on helicopter dynamics has been reported by Nagabhushanam et al. [3]; a review of various systems applied to structural dynamics has been given by Noor and Anderson [4]. A layman's guide to computer algebra and its applications was given by Pavelle et al. [5].
There are quite a number of systems available, notably FORMAC, MACSYMA, REDUCE, SCRATCHPAD and most recently SMP. Some of these were reviewed by Jensen and Niordsen [6]. The present authors have first hand experience only of REDUCE and our limited knowledge of symbolic computing indicates that this system is the most widely used in engineering applications. Some recent papers have appeared in which problems very similar to our own have been treated. Kiessling [7, 8] discussed REDUCE in the context of wind turbine problems, but unfortunately both papers appeared in rather obscure publications. Reference [8] is a yery thorough description of the use of symbolic computing in this type of work. Koppens [9] has considered some rotorcraft problems, less complex than those considered by Kiessling or in the present paper, but he retained the non-linear terms which, of course, is no hardship when using the computer to perform the algebra.
The growing interest in symbolic computing is aptly demonstrated by 'the advent of a new technical journal: Journal of Symbolic Computation. In a paper published there Fitch [1o] has summed up its position in the engineering field: “"(We) ... undertook a programme for the education of ... companies in both Sweden and Britain in computer algebra, and we discovered that once these techniques had been shown to them they immediately recognized that they would be able to ask questions that they had been suppressing as incapable of solution. Many of these companies are now acquiring algebra systems. At present it is unfortunately still the case that much algebra is being done by hand'". That is the view of a computer scientist. We hope that the present paper will reinforce his view and encourage many other engineers to use this powerful aid to analysis.
# 2. WIND TURBINE DYNAMICS
The main problem in producing a mathematical model of a wind turbine is the fact that the rotor exercises gross angular movements with respect to the support structure. In.standard structural analysis techniques, such as finite element packages, it is assumed that the structure has a mean position about which deflections, albeit large deflections, occur. This is evidently not the case for a wind turbine. The analyst is therefore forced to develop the equations of motion from first principles. The goal of the analysis is to be able to predict the stability of the system and to be able to calculate the forced response.
# 3. DERIVATION OF THE EQUATIONS OF MOTION
3.1. DEFINITION OF THE CO-ORDINATE SYSTEMS
It is common practice when faced with the derivation of the equations of motion of a complex dynamic system to use the well-known approach due to Lagrange. This requires the calculation of an expression for the various forms of energy in the system, the kinetic energy $\pmb{T,}$ the strain energy $\pmb{U}$ and the dissipation energy $U_{D}$ . After having calculated these quantities the Lagrangian may be derived and equated with an expression for the generalized force $Q.$
For a system such as we are considering the computation of the strain and dissipation energies is relatively straightforward compared with the kinetic energy and generalized force.
For the wind turbine problems we considered standard methods were used to set up the description of an arbitrary point on one of the blades. A series of transformations (rotations and translations) were calculated that allowed the position vector of an arbitrary point on the blades to be expressed in inertially fixed axes. These transformations are listed below and sketched in Figure 1.
The inertially fixed co-ordinate system is defined such that it is aligned along the unperturbed tower. The remaining co-ordinate systems are defined as follows: $\begin{array}{r}{T_{0}(x_{0},y_{0},}\end{array}$ $\scriptstyle z_{0})$ is the inertially fixed system; $T_{1}(x_{1},\,y_{1},\,z_{1})$ is obtained by translation of $\boldsymbol{T_{0}}$ by the tower head displacement vector, $\mathfrak{n}_{1};\;T_{2}(x_{2},\;y_{2},\;z_{2})$ is obtained by torsional rotation, $\theta_{y},$ about $y_{1};\,T_{3}(x_{3},\,y_{3},\,z_{3})$ is obtained by “"nodding" rotation, $\pmb{\theta_{x}}$ , about $x_{2};\,T_{4}(x_{4},\,y_{4},\,z_{4})$ is obtained by lateral rotation, $\pmb{\theta}_{\pmb{z}}$ , about $z_{3};\,I_{5}(x_{5},\,y_{5},\,z_{5})$ is obtained by translation of $\pmb{T_{4}}$ by the overhang vector $\mathbf{n_{2}}$ $T_{6}(x_{6},y_{6},z_{6})$ is obtained by azimuthal rotation, $(\psi+\phi_{z})$ ,about $z_{5};\,I_{7}(x_{7},\,y_{7},\,z_{7})$ is obtained by teeter rotation, $\pmb{\gamma},$ about $x_{6}$
当面临复杂动力系统运动方程的推导时,通常采用著名的拉格朗日方法。这需要计算系统中各种形式的能量表达式,即动能 $\pmb{T}$、应变能 $\pmb{U}$ 和耗散能 $U_{D}$。在计算出这些量之后,可以推导出拉格朗日量,并将其与广义力 $Q$ 的表达式相等。
对于我们所考虑的系统,与动能和广义力相比,应变能和耗散能的计算相对简单。
对于我们考虑的风电机组问题采用标准方法来建立叶片上任意点的描述。计算了一系列变换旋转和平移使得叶片上任意点的位矢能够在惯性固定坐标系中表示。这些变换列举如下并在图1中示意。
惯性固定坐标系定义为沿未受扰动的塔架对齐。其余坐标系定义如下:$T_{0}(x_{0},y_{0},z_{0})$ 是惯性固定坐标系;$T_{1}(x_{1},\,y_{1},\,z_{1})$ 是通过将 $\boldsymbol{T_{0}}$ 沿塔顶位移矢量 $\mathfrak{n}_{1}$ 平移得到;$T_{2}(x_{2},\;y_{2},\;z_{2})$ 是通过绕 $y_{1}$ 轴进行扭转旋转 $\theta_{y}$ 得到;$T_{3}(x_{3},\,y_{3},\,z_{3})$ 是通过绕 $x_{2}$ 轴进行“点头”旋转 $\pmb{\theta_{x}}$ 得到;$T_{4}(x_{4},\,y_{4},\,z_{4})$ 是通过绕 $z_{3}$ 轴进行侧向旋转 $\pmb{\theta}_{\pmb{z}}$ 得到;$I_{5}(x_{5},\,y_{5},\,z_{5})$ 是通过将 $\pmb{T_{4}}$ 沿悬伸矢量 $\mathbf{n_{2}}$ 平移得到;$T_{6}(x_{6},y_{6},z_{6})$ 是通过绕 $z_{5}$ 轴进行方位角旋转 $(\psi+\phi_{z})$ 得到;$I_{7}(x_{7},\,y_{7},\,z_{7})$ 是通过绕 $x_{6}$ 轴进行摆动旋转 $\pmb{\gamma}$ 得到。
![](images/3cc981e37406bc5eafd73cf77240b5fcb6041716321a2bd522f65b4209f5b9d1.jpg)
Figure 1. (a) System degrees of freedom; $\pmb{\mathrm{\vec{z}_{1}}}$ ,tower head displacement; $\theta_{y}$ , tower torsion; $\theta_{\mathbf{x}}$ ,tower nodding; $\pmb{\theta}_{\pmb{z}}$ , tower lateral rotation; $\mathfrak{n}_{2}$ nacelle overhang; $\psi+\phi_{z}$ , azimuthal rotation; $\pmb{\gamma},$ rotor teeter; $\mathbb{W}_{1}$ ,rotor flapwise bending; $w_{2}$ rotor edgewise bending; (b) co-ordinate systems.
The azimuthal rotation, $(\psi+\phi_{z})$ is the sum of the steady state and power train perturbation components.
The rotational matrices defining the co-ordinate transformations are standard and may be easily determined from the definition of the co-ordinate systems given above.
Generalized co-ordinates were used that corresponded to the following physical degrees of freedom; rotor teeter; rotor azimuth and power train rotation; flapwise and edgewise displacement including arbitrarily large initial deflections; tower torsion; tower fore-aft and side-to-side displacement.
方位角旋转 $(\psi+\phi_{z})$ 是稳态分量和传动链扰动分量之和。
定义坐标变换的旋转矩阵是标准的,并且可以从上面给出的坐标系定义中很容易地确定。
采用了与以下物理自由度相对应的广义坐标:风轮跷跷板运动;风轮方位角和传动链旋转;包括任意大初始变形的挥舞和摆振位移;塔架扭转;塔架前后和侧向位移。
# 3.2. DERIVATION OF KINETIC ENERGY
The kinetic energy of the wind turbine structure may be written as
$$
T\!=\!\frac{1}{2}\int_{{\mathrm{rotor}}}\mathbf{V}_{6}^{\mathsf{T}}M_{R}\mathbf{V}_{6}\,\mathrm{d}r\!+\!\frac{1}{2}\!\!\int_{{\mathrm{tower}}}\mathbf{V}_{0_{T}}^{\mathsf{T}}[M_{T}]\mathbf{V}_{0_{T}}\,\mathrm{d}h,
$$
where $\mathbf{v}_{\mathfrak{o}_{\tau}}$ and $\mathbf{y}_{6}$ are the absolute velocity vectors of a point on the tower and on the rotor, expressed in $\scriptstyle{T_{0}}$ and $\boldsymbol{T_{6}}$ co-ordinate systems respectively. The tower distributed mass and inertia properties are contained in the matrix $M_{T}$
After considerable manipulation the velocity of a point on the rotor can be expressed as
$$
\begin{array}{l}{{{\bf{V}}_{6}\!=\![\,T_{\bar{\nu}}]^{\top}[\,T_{\theta_{c}}]^{\top}[\,T_{\theta_{s}}]^{\top}[\,T_{\theta_{\nu}}]^{\top}\bar{\bf{n}}_{1}\!+\![\,T_{\bar{\nu}}]^{\top}[\,T_{\theta_{c}}]^{\top}[\,T_{\theta_{s}}]^{\top}(\dot{\bf{0}}_{y}\times([\,T_{\nu}]{\bf{R}}_{\bar{\nu}}\!+\![\,T_{\theta_{x}}][\,T_{\theta_{z}}]{\bf{n}}_{2}))}}\\ {{\mathrm{~}\qquad+[\,T_{\bar{\nu}}\,]^{\top}[\,T_{\theta_{z}}]^{\top}(\dot{\bf{0}}_{x}\times([\,T_{\gamma}]{\bf{R}}_{\bar{\nu}}\!+\![\,T_{\theta_{c}}]{\bf{n}}_{2}))\!+\![\,T_{\bar{\nu}}]^{\top}[\,\dot{\bf{0}}_{z}\times([\,T_{\gamma}]{\bf{R}}_{\bar{\nu}}\!+\!{\bf{n}}_{2}))}}\\ {{\mathrm{~}\qquad+\dot{\bar{\Psi}}\times[\,T_{\gamma}]{\bf{R}}_{7}\!+\![\,T_{\gamma}](\dot{\bf{R}}_{7}\!+\!\dot{\gamma}\times{\bf{R}}_{7}),}}\end{array}
$$
where the rotational matrices $[T]$ are standard and may easily be determined from the definition of the co-ordinate systems.
After the rotor absolute velocity has been calculated in this way, the kinetic energy is derived by substitution of the generalized co-ordinates into equation (1). The resultant expression is not included in this text since it is extremely cumbersome and of no great physical significance.
# 3.3. DERIVATION OF THE GENERALIZED FORCE
The generalized force corresponding to the arbitrary generalized co-ordinate $\pmb q$ is given by
$$
Q=\int\mathbf{F}_{0}\cdot{\frac{\partial\mathbf{R}_{0}}{\partial q}}\,\mathrm{d}r,
$$
where $\mathbf{F_{0}}$ is the aerodynamic force vector acting at a point on the rotor defined by the displacement vector $\mathbf{R_{0}}$
The rotor-fixed force vector $\mathbf{F}_{7}$ is more convenient to work with and is derived from consideration of the local aerodynamic lift and drag loads resolved into the blade edgewise and flapwise directions and is given by
$$
\begin{array}{r}{\mathrm{d}\mathbf{F}_{7}\!=\!\left(\!\frac{1}{2}\rho U c(C_{L}U_{P_{7}}\!-\!C_{D}U_{T_{7}})\:\mathrm{d}r\!\right),}\\ {\mathrm{d}\mathbf{F}_{7}\!=\!\left(\!\frac{1}{2}\rho U c(C_{L}U_{T_{7}}\!+\!C_{D}U_{P_{7}})\:\mathrm{d}r\!\right),}\end{array}
$$
where $U_{T_{7}}$ and $U_{P_{7}}$ are the edgewise and flapwise components of the relative wind velocity vector $\mathbf{V}_{\mathtt{R e I}_{7}}$ expressed in $\scriptstyle{T_{7}}$ rotor-fixed axes, and $\bar{U^{=}}(U_{T_{7}}^{2}\!+U_{P_{7}}^{2})^{1/2}$
The relative wind velocity vector defines the velocity of the wind with respect to the structural velocity of the rotor and is computed from
$$
\Psi_{\mathrm{Rel}_{7}}\!=\![\,T_{\gamma}][\,T_{\bar{\psi}}][\,T_{\theta_{\varepsilon}}][\,T_{\theta_{x}}][\,T_{\theta_{y}}](\Psi_{w_{0}}\!-\!\dot{\bf R}_{0}),
$$
where $\mathbf{v}_{\mathsf{w}_{0}}$ is the absolute velocity of the wind expressed in inertial axes and the $\dot{\mathbf{R}}_{0}$ is the structural velocity.
The aerodynamic force vector $\mathbf{F_{0}}$ is derived by transforming $\mathbf{F}_{7}$ Thegeneralizedforce may then be computed from equation (2) after having first substituted into it the generalized co-ordinates.·
# 3.4. EQUATIONS FOR THREE-BLADED MACHINES
It is interesting to consider the slightly different but strongly related problem of a three-bladed rotor. The discussion above described the derivation of the equations of motion of a two-bladed rotor starting with an arbitrary point on the blade. The teeter freedom requires the two-bladed rotor to be treated as a single entity rather than as two independent blades. In the basic steps involved an arbitrary point was considered and it is only the specification of the mode shapes themselves that narrows the analysis to a two-bladed case. One may, therefore, take the kinetic energy derivation of a blade on a flexible tower as a starting point for the analysis of a three-bladed machine. It was mentioned above that two-bladed rotors pose considerable numerical problems in stability calculations. A three-bladed rotor is much simpler since the co-ordinates may be transformed to remove the periodic terms. Dugundji and Wendell [11] have given a good review of this aspect of wind turbine dynamics. In our case we were able to make uise of the single-blade analysis that we had derived for the two-bladed work. Like the derivation of the equations of motion themselves, the transformation from physical co-ordinates into multi-blade co-ordinates can be very simply specified and was therefore suitable for symboliccomputing.
To illustrate this transformation consider a state vector that describes a three-bladed rotor supported by a flexible tower. The blades are modelled in terms of their first flapping and lead-lag modes. $\pmb q_{1,i}$ denotes the generalized co-ordinate of the first flapping mode of bladei and $\pmb{p}_{1,i}$ the first lead-lag co-ordinate.
$\phi_{z}$ is the power train angular deflection, $\pmb{n},\ \pmb{s},\ \pmb{\theta}_{\pmb{y}}$ are the fore-aft, side-to-side and torsional generalized co-ordinates of the tower. The state vector of the system is $\pmb{\mathrm{x}}$ where
$$
\mathbf{x}^{\mathsf{T}}\!=\![q_{1,1}\quad q_{1,2}\quad q_{1,3}\quad p_{1,1}\quad p_{1,2}\quad p_{1,3}\quad\phi_{z}\quad n\quad s\quad\theta_{y}].
$$
The system can be transformed into multi-blade co-ordinates by making standard substitutions and defining a new state vector ${\bf x^{i}}$
$$
\mathbf{x}^{1\!\mathrm{T}}\!=\![q_{0}\quad q_{\!}\quad q_{\!}\quad p_{0}\quad p_{s}\quad p_{c}\quad\phi_{z}\quad n\quad s\quad\theta_{y}].
$$
If the untransformed equation of motion was
$$
[M]{\ddot{\mathbf{x}}}+[C]{\dot{\mathbf{x}}}+[K]\mathbf{x}=\mathbf{0},
$$
where $[M],[C]$ and $[\kappa]$ are periodic, then the new equation of motion is
$$
[M^{1}]{\ddot{\bf x}}^{1}\!+\![C^{1}]{\dot{\bf x}}^{1}\!+\![K^{1}]{\bf x}^{1}\!=\!{\bf0},
$$
where
$$
\begin{array}{r}{[M^{1}]\!=\![P]^{\top}\![M][P],\qquad[C^{1}]\!=\![P]^{\top}\!(2[M][\dot{P}]\!+\![C][P]),}\\ {[K^{1}]\!=\![P]^{\top}\!([M][\ddot{P}]\!+\![C][\dot{P}]\!+\![K][P]),~~~~~~~~~~~~}\end{array}
$$
which contain only constant terms due to the resulting trigonometric summations. Here
$$
[P]\!=\!\!{\left[\begin{array}{l l l}{[Q]}&{0}&{0}\\ {0}&{[Q]}&{0}\\ {0}&{0}&{[I]}\end{array}\right]},\qquad[Q]\!=\!{\left[\begin{array}{l l l}{1}&{\sin\psi_{1}}&{\cos\psi_{1}}\\ {1}&{\sin\psi_{2}}&{\cos\psi_{2}}\\ {1}&{\sin\psi_{3}}&{\cos\psi_{3}}\end{array}\right]},
$$
Where $[I]$ is the $\pmb{4\times4}$ identity matrix and the variable $\psi_{i}$ is the azimuth of the ith blade.
3.5. SYSTEMATIC DERIVATION OF THE LAGRANGIAN AND GENERALIZED FORCE
Equation (1) appears fairly innocent, but reference to equation (2) which shows the algebraicformof $\mathbf{v_{6}}$ demonstrates that calculation of ${\cal T}_{R}$ is formidable indeed. The complexity and tedium of the algebra involved becomes still worse when the derivation of the equations of motion is considered and the evaluation of the Lagrangian is attempted. The approach adopted to derive the expression for $T_{R}$ is fairly standard and follows' very closely that described by Ottens and Zwaan [12]. One attempts to minimize the algebraic manipulation by avoiding the calculation of the position vector and then taking its time derivative to calculate the velocity which would be the most direct method. The complexity is reduced by calculating the inertial velocity in an intermediate set of axes, the result being $\mathbf{y}_{6}$ . Nevertheless, the algebra involved is extremely tedious and, in the view of the present authors, approaches the limits of human endurance or at least reliability. This sort of manipulation can, of course, be performed, but examination of the resulting expressions, which must be checked, suggested that short of performing a wholly independent analysis and comparing the results there was no possible, reliable method of checking.
The derivation of the kinetic energy and the steps involved in the derivation of the equations of motion are very well defined. It is possible to outline a general scheme for calculation of the Lagrangian in a step by step method: (1) specify the position vector $(\mathbf{\mathbf{R}}_{7})$ of an arbitrary point of the blade in blade axes; (2) specify the series of transformations that are required to express the position vector in inertial axes, ${\bf R}_{0}\!=\![T]{\bf R}_{7}$ (3) take the derivativeof $\scriptstyle\mathbf{R_{0}}$ with respect to time in order to calculate the inertial velocity, $\dot{\mathbf{R}}_{0}$ (4) calculate the kinetic energy of the point mass and integrate over the rotor, $\int m\dot{\bf r}_{i}^{2}\,{\bf d}r=T_{R}$ (5) form the Lagrangian for each of the generalized co-ordinates.
By following the five steps specified above it is possible to calculate the equation of motion for each of the generalized co-ordinates. Thus, after having completed step 4, step 5 may be repeated for each of the co-ordinates so that the complete system equations are produced. In the list above a very inelegant method has been adopted to calculate the velocity--the position vector is simply differentiated with respect to time. This approach makes the algebra still more tedious but makes the specification of the steps considerably simpler. Faced with the problem of checking the derivation of the equations of motion, and realizing the mechanical way in which they were derived, suggested to us that our problem was ideally suited to solution by means of an algebraic computing system.
As with the Lagrangian, the derivation of the generalized force vector appears to be a straightforward mathematical procedure. However, the algebraic complexity of the derivation is again a great obstacle to a manual solution. The problem is therefore well-suited to the application of symbolic computing techniques. It is possible to define a series of mathematical steps for calculation of the generalized force: (1) specify the wind velocity vector $\mathbf{V}_{\ast_{0}}$ in inertial axes and subtract the inertial rotor velocity, $\dot{\mathbf{R_{0}}}\mathbf{;}$ (2) compute the relative wind velocity in blade axes by means of the rotational transformations, $v_{\mathtt{R e l},}\!=\![T]^{\mathtt{T}}$ $(\bar{\mathbf{V}}_{\mathbf{w}_{0}}\!-\!\dot{\mathbf{R}}_{0})$ (see equation (5));(3) compute the aerodynamic force vector in blade axes by meansof $\mathbf{v}_{\mathtt{p e l},}$ and the lift and drag coefficients; (4) express the aerodynamic force vector in inertial axes, $\mathbf{F}_{0}\!=\![\,T]\mathbf{F}_{7};$ (5) take the partial derivative of $\mathbf{R_{0}}$ with respect to the generalized co-ordinate to give the vector, $\partial\mathbb{R}_{0}/\partial q$ ; (6) form the dot product of $\mathbf{F_{0}}$ with $\partial\mathbb{R}_{0}/\partial q$ and integrate over the rotor to give the generalized force.
The algebra involved in the derivation of the generalized force vector is clearly complex and yet the mathematical procedure remains relatively straightforward. It is true to say, however, that the procedure requires a greater degree of careful consideration at each step compared to the easily automated derivation of the Lagrangian.
方程(1)看起来相当简单,但参考显示 $\mathbf{v_{6}}$ 代数形式的方程(2)表明 ${\cal T}_{R}$ 的计算确实非常艰巨。在考虑运动方程的推导以及尝试评估拉格朗日量时,所涉及代数的复杂性和繁琐性变得更糟。采用的推导 $T_{R}$ 表达式的方法相当标准,并与 Ottens 和 Zwaan [12] 所述的方法非常相似。人们试图通过避免计算位矢,然后对其进行时间求导来计算速度,从而最小化代数运算,尽管这将是最直接的方法。通过在中间坐标系中计算惯性速度来降低复杂性,结果为 $\mathbf{y}_{6}$。尽管如此,所涉及的代数运算极其繁琐,在本文作者看来,这接近了人类耐力或至少是可靠性的极限。这种操作当然可以执行,但对必须检查的所得表达式进行检查,表明除了进行完全独立的分析并比较结果之外,没有其他可能可靠的检查方法。
动能的推导以及运动方程推导中涉及的步骤定义得非常清楚。可以概述一种分步计算拉格朗日量的通用方案:(1) 在叶片坐标系中指定叶片任意点的位矢 $(\mathbf{\mathbf{R}}_{7})$(2) 指定将位矢表示为惯性坐标系所需的系列变换,${\bf R}_{0}\!=\![T]{\bf R}_{7}$(3) 对 $\scriptstyle\mathbf{R_{0}}$ 进行时间求导以计算惯性速度 $\dot{\mathbf{R}}_{0}$(4) 计算质点的动能并在风轮上积分,$\int m\dot{\bf r}_{i}^{2}\,{\bf d}r=T_{R}$(5) 为每个广义坐标形成拉格朗日量。
通过遵循上述五个步骤可以计算每个广义坐标的运动方程。因此在完成步骤4后可以对每个坐标重复步骤5以便生成完整的系统方程。在上述列表中采用了一种非常不优雅的方法来计算速度——即简单地对位矢进行时间求导。这种方法使代数运算更加繁琐但使步骤的规定大大简化。面对检查运动方程推导的问题并意识到它们的机械推导方式使我们认为我们的问题非常适合通过代数计算系统来解决。
与拉格朗日量一样,广义力矢量的推导似乎是一个直接的数学过程。然而,推导的代数复杂性再次成为手动解决的巨大障碍。因此,该问题非常适合应用符号计算技术。可以定义一系列计算广义力的数学步骤:(1) 在惯性坐标系中指定风速矢量 $\mathbf{V}_{\ast_{0}}$ 并减去惯性风轮速度 $\dot{\mathbf{R_{0}}}$(2) 通过旋转变换计算叶片坐标系中的相对风速,$v_{\mathtt{R e l},}\!=\![T]^{\mathtt{T}}$ $(\bar{\mathbf{V}}_{\mathbf{w}_{0}}\!-\!\dot{\mathbf{R}}_{0})$ (参见方程(5))(3) 通过 $\mathbf{v}_{\mathtt{p e l},}$ 以及升力系数和阻力系数计算叶片坐标系中的气动力矢量;(4) 将气动力矢量表示为惯性坐标系,$\mathbf{F}_{0}\!=\![\,T]\mathbf{F}_{7}$(5) 对 $\mathbf{R_{0}}$ 关于广义坐标求偏导,得到矢量 $\partial\mathbb{R}_{0}/\partial q$(6) 将 $\mathbf{F_{0}}$ 与 $\partial\mathbb{R}_{0}/\partial q$ 进行点积,并在风轮上积分以得到广义力。
广义力矢量推导中涉及的代数运算显然很复杂,但数学过程仍然相对直接。然而,可以说,与易于自动化的拉格朗日量推导相比,该过程在每个步骤都需要更仔细的考虑。
# 4.SOME EXAMPLES OF SYMBOLIC MANIPULATIONS
Symbolic computing comes into its own when tackling complicated and tedious prob. lems. The type of problem for which we have found it useful is described in section 3. It will come as no surprise to workers in the field of dynamics that the matrices involved in these problems are very large and, although it is vital that they are correct, it is difficult, if not impossible, to glean any physical understanding of the system they describe by scrutinizing them term by term.
The problems described above have been analyzed by using symbolic computing techniques but, rather than tackle a large problem in detail here, which will result in large and complex algebraic expressions, it is more illustrative to consider a very simple example which may be followed through from beginning to end.
![](images/a2cd117bf0a8fcb5957e734db7f3e18f6c3f3d4c0b0c9f1daba042d2d36f4a6e.jpg)
Figure 2. Three degree of freedom model. $\psi=\pmb{\Omega}t.$
Consider the wind turbine model shown in Figure 2. This is a simple three degree of freedom model first described by Kaza et al. [13]. It bears little physical resemblance to a real wind turbine, in fact its dynamic characteristics correspond rather more closely to a helicopter; however, it does do well as an illustrative model. The degrees of freedom contained in the model are tower head lateral motion and blade lead-lag motion. The derivation of the equations of motion of this system follows exactly the same steps as for the larger system described in section 3. The equivalent REDUCE steps, together with the inertial parts of the equations themselves, are given in Table 1.
TABLE 1 A simple example of REDUCE
$\%$ The following REDUCE program calculates the Lagrangians of
$\%$ the simple 3 d.o.f system described in Figure 2.
$\%$
MATRIX MPSI(3,3),VTOW(3,1),VR1(3,1),VR2(3,1),VR3(3,1),VR3DOT(3,1)S $\%$
$\%$ Algebraic procedure for differentiation w.r.t. time:
$\%$
PROCEDURE DIFFRT(A)S
BEGIN
RETURN QDOT\*DF(A,Q)
+PSIDOT\*DF(A,PSI)S
ENDS
$\%$
$\%$ Algebraic procedure for determination of a Lagrangian:
$\%$
PROCEDURE LAG(B,C)S
BEGIN
$\begin{array}{r}{\mathsf{D D}\!:=\!\mathsf{D F}(\mathsf{T F},\!\mathbf{B})\mathbb{S}}\\ {\mathsf{A A}\!:=\!\mathsf{D F}(\mathsf{T F},\!\mathbf{C})\mathbb{S}}\end{array}$
RETURN OMEGA\*DF(AA,PSI0)
+QDOT\*DF(AA,Q)
+QDOTDOT\*DF(AA,QDOT)
+Z1DOT\*DF(AA,Z1)
+Z1DOTDOT\*DF(AA,Z1DOT)
+Z2DOT\*DF(AA,Z2)
+Z2DOTDOT\*DF(AA,ZDOT)
-DDS
ENDS
$\%$
$\%$ MPSI is the azimuthal rotation matrix:
$\%$
MPSI := MAT(COS(PSI),-SIN(PSI),0),(SIN(PSI),COS(PSI),0),(0,0,1))S
$\%$
$\%$ VTOw is the tower head displacement defined in inertial axes:
$\%$
$\mathbf{VIOW}{:=}\,\mathbf{MAT}((\mathbf{Q}),(\varnothing),(\varnothing))\mathbb{S}$
$\%$
$\%$ VR1 is the position vector of an arbitrary point on the blade in blade axes: $\%$
$\mathbf{VR1}\!:=\mathbf{MAT}((\boldsymbol{\mathfrak{G}}),(\mathbf{RR}),(\boldsymbol{\mathfrak{G}}))\mathbf{S}$
VR2:= MPSI\*VR1S
$\begin{array}{r}{\mathbf{V}\mathbf{R}3:=\mathbf{V}\mathbf{R}2\mathbf{+}\mathbf{V}\mathbf{T}0\mathbf{W}\S}\end{array}$
CLEAR MPSI,VTOW,VR1, VR2S
$\%$
$\%$ VR3DOT is the general blade velocity vector defined in inertial axes:
$\%$
VR3DOT $:=$ MAT(DIFFRT(VR3(1,1)),(DIFFRT(VR3(2,1)),(DIFFRT(VR3(3,1)s CLEAR VR3S
$\%$
$\%$ VSQ is the velocity squared expression for an arbitrary blade station:
$\%$
VSQ := VR3D0T(1,1)\*VR3D0T(1,1) $^{1+}$ VR3DOT(2,1)\*VR3DOT(2,1)+
VR3DOT(3,1)\*VR3D0T(3,1)S
LET SIN(PSI) $\lvert=$ SINPSIS
LET COS(PS $\scriptstyle\mathtt{I}\mathtt{J}=\mathtt{I}$ COSPSIS
LET $\mathrm{COSPSI^{**}2}\!\!=\!1\!-\!\mathrm{SINPSI^{**}2S}$
$\mathbf{VSQ:=VSQS}$
TABLE 1 (cont.)
LET SINPSI\*\*2=1-COSPSI\*\*2S
$\mathbf{VSQ:=VSQ5}$
$\%$ Z1 and ${\bf z}{\bf2}$ are the perturbation co-ordinates for blades I and 2:
$\%$
PSIDC $)\mathbf{T}{:=}\,\mathbf{OMEGA}$ +Z1DOTS
LET $\mathbf{SIN}(\mathbf{ZI})\!=\!\mathbf{Z1}\mathfrak{S}$
LET $\mathsf{C O S}(\mathbf{Z1})\!=\!1\!-\!z1^{**}\!2/2\mathfrak{I}$
SINPSI := SIN(PSI0)\*COS(Z1)+COS(PSI0)\*SIN(Z1)\$
$\mathbf{COSPSI:=COS(PSI/\beta)^{\ast}C O S(Z1)\mathbf{-SIN(PSI/\beta)^{\ast}S I N(Z1)\mathfrak{H}}}$
$\%$
$\%$ VSQ1 and VSQ2 are the velocity squared expressions for blades 1 and 2: $\%$
$\mathbf{VSQ1:=VSQ95}$
$\mathbf{VSO2:=SUB}$ (Z1=Z2,Z1DOT $\ddots$ Z2DOT,SIN(PSI0) $=$ -SIN(PSI0,COS(PSI0)= $-\mathbf{COS}(\mathbf{PSI\emptyset}),\mathbf{vSQ1})\mathbb{S}$
VSQTOT $:=$ VSQ1+VSQ2\$
CLEAR VSQ,VSQ1,VSQ2\$
$\%$
$\%$ TT is the total kinetic energy of the rotor:
$\%$
TT:= M\*VSQTOT/2\$
$\%$
$\%$ Definition of blade inertia IB and first mass moment integral SB:
$\%$
LET M\*RR\*RR $\cdot\overline{{-}}$ IB\$
$\mathbf{T}\mathbf{T}:=\mathbf{T}\mathbf{T}\mathbf{\mathbb{S}}$
LET $\Delta A^{*}\mathbb{R}\mathbb{R}{=}\mathbb{S}\mathbb{B}\mathbb{S}$
$\mathbf{T}:=\mathbf{T}\mathbf{T}\mathbf{S}$
CLEAR Q,QDOT,Z1,Z1DOT,Z2,Z2DOT\$
$\%$
$\%$ The Lagrangians are now computed:
$\%$
$\mathbf{LAGQ}:=\mathbf{LAG}(0,\mathbf{QDOT})\mathbf{\hat{S}}$
$\mathbf{LAGZ1\!:=\!LAG(Z1,\!Z1DOT)\!\!\!\S}$
$\mathbf{LAGZ2}:=\mathbf{LAG}(\mathbf{Z}2,\mathbf{Z2DOT})\hat{\mathbf{S}}$
$\%$
$\%$ A weight level is defined for linearization of the Lagrangians.
$\%$ The linearized expressions are then output;
$\%$
WTLEVEL 1\$
WEIGHT $\mathbf{Q}\mathbf{=}1,\mathbf{Q}\mathbf{D}\mathbf{O}\mathbf{T}\mathbf{=}1,\mathbf{Q}\mathbf{D}\mathbf{O}\mathbf{T}\mathbf{D}\mathbf{O}\mathbf{T}\mathbf{=}1,\mathbf{Z}1\mathbf{=}1,\mathbf{Z}1\mathbf{D}\mathbf{O}\mathbf{T}\mathbf{=}1,\mathbf{Z}1\mathbf{D}\mathbf{O}\mathbf{T}\mathbf{D}\mathbf{O}\mathbf{T}\mathbf{=}$
$1,\!Z\!2\!=\!1,\!Z\!2\mathrm{DOT}\!\!=\!1,\!Z\!2\mathrm{DOTDOT}\!\!=\!1\mathfrak{F}$
$\mathbf{LAGQ}:=\mathbf{LAGQ};$
$\mathbf{LAGZ1}:=\mathbf{LAGZ1}:$
$\mathbf{LAGZ2}:=\mathbf{LAGZ2}$
ENDS
$\%$ Lagrangian for 'q generalized coordinate;
$\mathrm{LAGQ}:=\mathrm{COS}(\mathrm{PSI}\emptyset)^{\ast}\mathrm{Z2DOTDOT}^{\ast}\mathrm{SB}-\mathrm{COS}(\mathrm{PSI}\emptyset)^{\ast}\mathrm{Z2}^{\ast}\mathrm{OMEGA}^{\ast}\mathrm{SB}$
-COS(PSI0)
$\mathrm{^{*}Z1D O T D O T^{*}S B+C O S(P S I/6)^{*}Z1^{*}O M E G A^{*}S B-2^{*}S I N(P S I/6)^{*}}$
$\mathrm{Z2DOT^{+}O M E G A^{*}S B+2^{*}S I N(P S I\beta)^{*}Z1D O T^{*}O M E G A^{*}S B+2^{*}M^{*}Q D O T D O T^{-}}$ $\%$ Lagrangian for Z1' generalized coordinate;
$\mathrm{LAGZ1}:=-\mathrm{COS}(\mathrm{PSI}\emptyset)^{\ast}\mathrm{QDOTDOT^{\ast}S B}+\mathrm{Z1DOTDOT^{\ast}I B}$
$\%$ Lagrangianfor $\mathbf{\mathcal{Z}}2^{\bullet}$ generalized coordinate: LAGZ2 := COS(PSI0)\*QDOTDOT\*SB+ Z2DOTDOT\*1B
The steps described in the table are fairly straightforward for any reader familiar with FORTRAN and dynamics. Some general comments may, however, be made. REDUCE has the capability of implementing "Procedures" that are very much akin to FORTRAN subroutines. In the example here we have used two procedures, one for differentiation with respect to time and one to form the Lagrangian as defined in equation (1). Note that REDUCE can only perform partial differentiation. The derivation of the equations of motion involves considerable manipulation of matrices. The simple example cited here shows how REDUCE can perform these operations. It is often necessary, even when manipulation is being performed automatically, to linearize expressions. This may simply be done by REDUCE by using the “"WEIGHT' statement. This statement is used at the end of the example to assign weights to the individual variables and, given a general weight level of 1, all quadratic and higher order terms are discarded. The resulting Lagrangians for each of the three degrees of freedom are then calculated and printed at the end of the example.
The WEIGHT command may also be used to reduce the size of the kinetic energy expression. When considering a particular generalized co-ordinate many of the terms are redundant since they disappear when the expression is differentiated. Judicious use of the WEIGHT and WTLEVEL statements avoids the computation of the redundant terms which brings about significant reductions in the size of the expressions used. Failure to carry out this procedure does, in our experience, lead to unacceptably large expressions which brings the calculation to a halt when storage is exhausted. This problem of mid-calculation “expression swell is very common and has been discussed by Nagabhushanam et al. [3].
One further, and very useful, facility that REDUCE possesses is the ability to write expressions to files in FORTRAN format thus enabling the computer to both derive the equation of motion and transcribe them into a FORTRAN program for numerical solution. This final step helps to eliminate further possible sources of careless mistakes, since a manual transcription process is a very likely place for such problems to arise.
# 5. SOME COMMENTS ABOUT THE USE OF REDUCE
The automatic derivation of the equations of motion of large dynamical systems brings with it two main advantages: increased confidence in the results and increased scope in the scale of problems that may be tackled. Some limitations are still imposed either by financial constraints or limitations in computer storage capacity, but, in principle, the size and complexity of a mathematical model need not be limited by the size and complexity of the algebraic expressions which it produces.
The analyses described here exercized only a small part of the REDUCE system. We have used matrix manipulation and differential calculus. The system contains a great many more analytical facilities particularly in the field of high energy physics. For the engineering user the integration and factorization facilities are probably the most relevant. It is interesting to note that the integrator does not simply run through a “"look-up table" of known integrals but tackles each one on a rigorous and systematic basis. There is still some considerable work to be done on this aspect of the system; for example, integrals containing square roots present great difficulties and REDUCE often fails to produce an answer. In view of the indirect nature of integration this is not surprising. Much work is evidently underway on this subject at present.
We have found REDUCE very useful. We have also found it to be difficult to implement and sometimes very frustrating to use. It would be erroneous to give the impression that problems such as those described in this paper have been solved without difficulty. The highly sophisticated nature of REDUCE which gives it its wide applicability and power brings with it much scope for strange errors and obtuse messages. However the new release, REDUCE 3·2 is evidently more robust and has overcome many of these problems.
We approached symbolic computing in a cautious frame of mind. We first attempted to do small problems with solutions that could also be obtained manually with reasonable ease. We then proceeded to use REDUCE to check a very long and complicated analysis as described above. Only when we had achieved term for term agreement with this calculation did we feel happy to use the technique alone. Although there are certain features of our present system (REDUCE 3·0 implemented on an IBM 3033 operating CMS running under VM/SP release 2) that do not work, we have not discovered any way of obtaining erroneous results. Discussions with other users have also failed to bring any such problems to light. Our careful approach was therefore unnecessary which is no surprise since REDUCE has been in use in high energy physics for many years.
# 6. RESULTS FROM THE WIND TURBINE CALCULATIONS
It is hoped that readers of this paper will be interested in both the application of symbolic computing to engineering problems and in the problems themselves. The original reason for embarking on the analyses described here was to ascertain the stability characteristics of wind turbine designs. It is well known that helicopters, which have some similarity with wind turbines, are prone to severe mechanical and aeroelastic instabilities. Wind turbines seem less susceptible to these problems, but it is nevertheless prudent to check for possible unstable regions within the operational envelope.
It was pointed out earlier that the simple three degree-of-freedom system used as an example above had dynamic characteristics which had more in common with helicopters than wind turbines. This comment is confirmed by consideration of the results presented in Figure 3. These were computed by using standard Floquet techniques together with some more economical methods; these calculations have been described by Quarton and Garrad [14]. They clearly show a region of instability. Investigation of the instability reveals it to be a “ground resonance" type. Such an instability requires the blades to be very soft in the lead-lag direction. The present generation of wind turbines tend to have blades that are very stiff in this direction and hence this type of instability cannot occur.
![](images/c8d820d96f90b393f0c6fe9aa78bd14a03824908d65c486af171f72eae40e405.jpg)
Figure 3. System frequencies of three degree of freedom model. -——, Rotor collective mode; -- -, rotor cyclicmode;·..,tower lateral mode.
![](images/c3a00899a6eb2ce0606a1aba87996943bfbfb45ce4b41d11362028aa0bdae146.jpg)
Figure 4. System WEG MS-2 frequencies.
![](images/f1ce29a078c5953544017d5a415330c0c1322ee43b8a0abd49e258cc8e23f18b.jpg)
Figure 5. Predicted and measured WEG MS-1 blade bending moments at $16\,\Pi/\mathfrak{s}$ windspeed.-, Predicted; -- -, measured. (a) Edgewise bending moment at $3{\cdot}28\,\mathrm{m}$ ; (b) fapwise bending moment at $\mathbf{3\cdot28}\,\mathbf{m}$
A more typical set of stability results is shown in Figure 4. These were obtained for the WEG MS-2 wind turbine (a $25\,\uppi$ diameter,3-bladed, $200\,\mathbf{kW}$ machine). This model contains the fundamental fapping and lagging modes of each blade, a rigid body drive train freedom, and the first fore-aft, side-to-side and torsional modes of the tower. The eigenvalues of the system were obtained from a constant coefficient set of equations which resulted from performing a Coleman multi-blade transformation on the set of periodic coefficient equations as described above in section 3.4. The equations of motion of the system were derived entirely by REDUCE which also wrote a considerable part of the FORTRAN program used to calculate the numerical results.
In addition to using REDUCE for the stability analysis of wind turbine systems it has also been used for the derivation of mathematical models for the calculation of forced response. The procedure for deriving the equations of motion for a forced response model is essentially the same as for stability analysis except it then becomes important to retain all steady forcing and defection terms as well as those which are proportional to the system generalized co-ordinates.
A forced response model for the WEG $_{20\,\mathfrak{m}}$ diameter,2-bladed, $250\,\mathbf{kW}$ MS-1wind turbine has been derived by using REDUCE and typical response predictions are presented in Figure 5. These predictions of blade bending moment have been compared with equivalent measured results and their satisfactory agreement provides validation of the mathematicalmodel.
# 7. CONCLUDING REMARKS
The purpose of this paper was to describe one area of engineering analysis in which symbolic computing has played a useful role. It is hoped that the examples described here together with the specific comments on the use of the REDUCE system will enable other analysts to consider the possibility of using symbolic computing to help tackle complicated algebraic problems.
# ACKNOWLEDGMENTS
The authors wish to thank the U.K. Department of Energy and the Directors of Taylor Woodrow Construction, both of whom have supported parts of this work, for their permission to publish this paper. We have received invaluable help over our problems with REDUCE from Dr Pat Pearce of Kingston Polytechnic and Professor John Fitch of the University of Bath. We are also very much in debt to Rob Warren of Taylor Woodrow Construction and Dr Larry Seward of the Rand Corporation, who have spent a lot of time and energy implementing the REDUCE system on Taylor Woodrow's IBM3033.
# REFERENCES
1. H.G. KAHRIMANIAN 1953 MS Thesis, Temple University, Philadelphia, Pennsylvania, Analytical differentiation by digital computer.
2. J. NOLAN 1953 SM Thesis, MIT Cambridge, Massachusetts. Analytical differentiation on a digital computer.
3. J. NAGABHUSHANAM, G. H. GAONKAR and T. S. R. REDDY 1981 Proceedings of the Seventh European Rotorcraft and Powered Lift Aircraft Forum, Garmisch-Partenkirchen,Federal Republic of Germany, Paper No 37. Automatic generation of equations for rotor body systems with dynamic inflow for a priori ordering schemes.
4. A. K. NoOR and C. M. ANDERsON 1979 Computers and Structures 10, 95-118. Computerised symbolic manipulation in structural mechanics progress and potential.
5. R. PAVELLE, M. ROTHSTEIN and J. FITCH 1981 Scientific American, 101-113. Computer algebra.
6. J. JENSEN and F. N1ORDsEN 1977 in Structural Mechanics Software Series (Editors N. Perrone and W. Pilkey) Charlottesville: University Press of Virginia. Symbolic and algebraic manipulation languages and their applications in mechanics.
7. F. KIEsSLING 1982 ICAS Proceedings, 13th Congress of the International Council of the Aeronautical Sciences/ AIAA Aircraft Systems and Technology Conference, Seattle, Washington. Computer-aided derivation of equations of motion for rotary-wing aeroelastic problems.
8. F. KIEssLING 1984 DFVLR Forschunsbevicht 84-10, 187 pages. Modelling of the overall aeroelastic system of a wind turbine aided by symbolic programming. (In German.)
9. W. P. KoPpENS 1983 Technische Hogeschool Delft, Luchtuaart en Ruimtevaarttechnik, Memorandum M477. Utilization of the symbolic manipulation system REDUCE 2 for deriving equations of motion.
10. J. FITCH 1985 Journal of Symbolic Computation 1, 211-227. Solving algebraic problems with REDUCE.
11. J. DUGUNDJ1 and J. H. WENDELL 1983 American Institute of Aeronautics and Astronautics Journal 21, 890-897. Some analysis methods for rotating systems with periodic coeffcients.
12. H. H. OTTENS and R. J. ZWAAN 1978 NLR TR 78115 L. Description of a method to calculate the aeroelastic stability of a two-bladed horizontal axis wind turbine.
13. K. R. V. KAZA, D. C. JANETZKE and T. L. SULLIVAN 1980 Journal of Energy 4, 162-169. MOSTAS computer code evaluation for dynamic analysis of two-bladed wind turbines. (Also NASA TM 79101.)
14. D. C. QUARTON and A. D. GARRAD 1984 Proceedings of the 6th British Wind Energy Association Conference, pp. 197-209. Some comments on the stability analysis of horizontal axis wind turbines.