vault backup: 2025-06-16 17:09:37

.obsidian/plugins/copilot/data.json

@@ -266,7 +266,7 @@
     },
     {
       "name": "Translate to Chinese",
-      "prompt": "<instruction>Translate the text below into Chinese:\n    1. Preserve the meaning and tone\n    2. Maintain appropriate cultural context\n    3. Keep formatting and structure\n    4. Blade翻译为叶片，flapwise翻译为挥舞，edgewise翻译为摆振，pitch angle翻译成变桨角度，twist angle翻译为扭角，rotor翻译为风轮，span翻译为展向，deflection翻译为变形，mode翻译为模态，normal mode翻译为简正模态，jacket 翻译为导管架，superelement翻译为超单元\n    Return only the translated text.</instruction>\n\n<text>{copilot-selection}</text>",
+      "prompt": "<instruction>Translate the text below into Chinese:\n    1. Preserve the meaning and tone\n    2. Maintain appropriate cultural context\n    3. Keep formatting and structure\n    4. Blade翻译为叶片，flapwise翻译为挥舞，edgewise翻译为摆振，pitch angle翻译成变桨角度，twist angle翻译为扭角，rotor翻译为风轮，turbine、wind turbine翻译为机组、风电机组，span翻译为展向，deflection翻译为变形，mode翻译为模态，normal mode翻译为简正模态，jacket 翻译为导管架，superelement翻译为超单元，shaft翻译为主轴\n    Return only the translated text.</instruction>\n\n<text>{copilot-selection}</text>",
       "showInContextMenu": true,
       "modelKey": "gemma3:12b|ollama"
     },

书籍/力学书籍/BladedTheoryManual/Model Linearisation/auto/Model Linearisation.md

@@ -1,7 +1,7 @@
 # Linear Analysis Background  
 
 The linear analysis calculations reduce the Bladed aeroelastic model to a linear system at each operating point requested by the user. The linear system of equations in state-space form is represented by  
-线性分析计算将Bladed气弹模型简化为用户请求的**每个工作点处的线性系统**。**以状态空间形式表示的线性方程组为**：
+线性分析计算将Bladed气弹模型简化为用户请求的**每个工作点处的线性系统**。**以状态空间形式表示的线性系统方程组为**：
 
 $$
 \begin{array}{r}{\underline{{\dot{\mathbf{x}}}}=\mathbf{A}\underline{{\mathbf{x}}}+\mathbf{B}\underline{{\mathbf{u}}}}\\ {\underline{{\mathbf{y}}}=\mathbf{C}\underline{{\mathbf{x}}}+\mathbf{D}\underline{{\mathbf{u}}}}\end{array}
@@ -27,7 +27,7 @@ It is important to note that in order to enable proper linearised wind turbine d
 · Physical effects that cannot be linearised shall be removed, for instance wind turbulence, stickslip, etc.  
 其中 $\mathbf{x}$ 是状态向量，代表系统状态；$\mathbf{u}$ 是系统输入向量；$\mathbf{y}$ 是系统输出向量。归一化向量 $\underline{{\mathbf{x}}}、\underline{{\mathbf{y}}}$ 和 $\underline{{\mathbf{u}}}$ 代表偏离平衡态的量。
 
-矩阵 $\mathbf{A}$、$\mathbf{B}$、$\mathbf{C}$ 和 $\mathbf{D}$ 代表这些向量之间的线性化关系。这代表对完整叶片（Blade）模型的一个简化，该模型使用一组完全非线性方程来计算状态导数和输出。
+矩阵 $\mathbf{A}$、$\mathbf{B}$、$\mathbf{C}$ 和 $\mathbf{D}$ 代表这些向量之间的线性化关系。这代表对完整（Bladed）模型的一个简化，该模型使用一组完全非线性方程来计算状态导数和输出。
 
 $$
 \begin{array}{r}{\dot{\mathbf{x}}=f(t,\mathbf{x},\mathbf{u})}\\ {\mathbf{y}=h(t,\mathbf{x},\mathbf{u}).}\end{array}
@@ -61,7 +61,7 @@ The elements of the matrices A, B, C and D can then be derived by performing a l
 
 在版本4.7及更早版本中，气动学状态不包含在模型线性化中。在版本4.8及更高版本的气动学公式中，用户可以选择是否包含气动学状态。
 
-为了进行线性分析，Bladed依次取每个工作点，并找到涡轮机的稳态条件（如同时域运行中的初始条件）。这意味着风轮没有加速，模态变形使得弹性载荷平衡外部载荷。这定义了$\mathbf{x_{0}},\,\mathbf{y_{0}}$和$\mathbf{u_{0}}$的值，即所有扰动围绕的主要平衡点。
+为了进行线性分析，Bladed依次取每个工作点，并找到风电机组的稳态条件（如同时域运行中的初始条件）。这意味着风轮没有加速，模态变形使得弹性载荷平衡外部载荷。这定义了$\mathbf{x_{0}},\,\mathbf{y_{0}}$和$\mathbf{u_{0}}$的值，即所有扰动围绕的主要平衡点。
 
 对于每个输入或状态，Bladed然后在平衡点两侧进行一系列幅度逐渐增加的扰动。人为地增加或减少状态或输入的数值，用这些修改后的数值求解系统，并记录状态导数和输出。扰动的数量和最大扰动幅度可以由用户定义。
 
@@ -258,7 +258,7 @@ $$
 ## Transforming the A,B,C,D matrices  
 
 We consider the linear model equations in the rotating frame of reference and define  
-我们考虑在旋转参考系中的线性模型方程，并定义：
+我们考虑在旋转参考系下的线性模型方程，并定义：
 $$
 \begin{align}
 
@@ -346,11 +346,11 @@ In their raw form, these eigenvector contributions represent the relative displa
 
 The phase of each contribution, $\phi_{i\,,}$ is determined by the argument of the corresponding complex eigenvector element, $v_{i},$ i.e.  
 
-每个耦合模态的未耦合模态贡献由其本征向量决定。如果耦合模态具有二阶状态（结构状态）的贡献，这些状态由状态向量中的两个状态表示，则使用位移状态来确定贡献。
+每个耦合模态的未耦合模态贡献由其特征向量决定。如果耦合模态具有二阶状态（结构状态）的贡献，这些状态由状态向量中的两个状态表示，则使用位移状态来确定贡献。
 
-以原始形式，这些本征向量贡献代表每个模态的相对位移，可用于构建耦合模态形状。然而，在坎贝尔图中，这些贡献已被归一化。这是通过修改本征向量矩阵来实现的，使得每一行和每一列都具有单位和。这会增加具有高质量和刚度的模态的百分比贡献，这些模态在位移中贡献很小，但在能量方面贡献很大。
+以原始形式，这些特征向量贡献代表每个模态的相对位移，可用于构建耦合模态形状。然而，在坎贝尔图中，这些贡献已被归一化。这是通过修改特征向量矩阵来实现的，使得每一行和每一列都具有单位和。这会增加具有高质量和刚度的模态的百分比贡献，这些模态在位移中贡献很小，但在能量方面贡献很大。
 
-每个贡献的相位，$\phi_{i\,,}$ 由相应复本征向量元素的参数决定，即 $v_{i},$ 即。
+每个贡献的相位，$\phi_{i\,,}$ 由相应复特征向量元素的参数决定，即 $v_{i},$ 即。
 $$
 \phi_{i}=\arctan\left(\frac{\operatorname{Im}\left(v_{i}\right)}{\operatorname{Re}\left(v_{i}\right)}\right).
 $$  
@@ -387,7 +387,7 @@ After the transformation and renaming of the individual blade modes the coupled
 
 Coupled mode name  1st uncoupled mode  2nd uncoupled mode  Phase angle $(\phi_{2}-\phi_{1})$   
 
-如果应用了 MBC 变换，则各个叶片模态将被转换成一组风轮模态。对于三叶风轮，通常存在集变模态collective、余弦循环模态cosine-cyclic和正弦循环sine-cyclic模态。所有叶片的 1 阶挥舞模态将被重命名为风轮 1 阶挥舞集变模态、风轮 1 阶挥舞正弦循环模态和风轮 1 阶挥舞余弦循环模态。如果叶片数量为偶数，则还会存在摆振模态。
+**如果应用了 MBC 变换，则各个叶片模态将被转换成一组风轮模态**。对于三叶风轮，通常存在集变模态collective、余弦循环模态cosine-cyclic和正弦循环sine-cyclic模态。所有叶片的 1 阶挥舞模态将被重命名为风轮 1 阶挥舞集变模态、风轮 1 阶挥舞正弦循环模态和风轮 1 阶挥舞余弦循环模态。如果叶片数量为偶数，则还会存在摆振模态。
 
 在转换和重命名各个叶片模态之后，命名耦合风轮模态。旋转模态的识别遵循下表中的逻辑。
 

书籍/力学书籍/BladedTheoryManual/Structural Dynamics/auto/Structural Dynamics.md

@@ -205,7 +205,8 @@ Last updated 26-11-2024
 
 A schematic of the Bladed calculation to evaluate the structural response in the time domain is shown in Figure 1. The numbered steps below map directly to the numbered steps in the diagram.  
 
-1. Modal displacements and velocities ("state values") are known at the start of the time step. 2. Applied loads are calculated based on the external loads and the state values. External loads applied at structural nodes are transformed into applied loads on the modes.  
+1. Modal displacements and velocities ("state values") are known at the start of the time step. 
+2. Applied loads are calculated based on the external loads and the state values. External loads applied at structural nodes are transformed into applied loads on the modes.  
 
 3. The structural dynamics equation $\mathbf{M}\ddot{\mathbf{q}}+\mathbf{C}\dot{\mathbf{q}}+\mathbf{K}\mathbf{q}=\mathbf{f}$ is solved in modal space to find the state accelerations q, where  
 
@@ -223,7 +224,7 @@ In most cases, applied loads depend partially on the nodal positions and velocit
 	-  $M、C、K$分别是系统模态矩阵，代表质量、阻尼和刚度；
 	-  $\mathbf{q}、\dot{\mathbf{q}}、\ddot{\mathbf{q}}$分别是状态位移、速度和加速度；
 	-  $f$是作用于每个状态的外部力模态向量。
-4. 积分器使用加速度来找到下一个时间步的状态值。
+4. 积分器==使用加速度==来找到下一个时间步的状态值。
 
 在大多数情况下，施加的载荷部分地取决于从状态值计算出的节点位置和速度，这很方便，因为状态值在每个时间步的开始时是已知的。
 

多体+耦合求解器/参考文献/Craig-Coupling of substructures for dynamic an/auto/Craig-Coupling of substructures for dynamic an.md

@@ -0,0 +1,592 @@
+# A00-24585  
+
+# AIAA-2000-1573  
+
+# COUPLING OF SUBSTRUCTURES FOR DYNAMIC ANALYSES: AN OVERVIEW  
+
+Roy R. Craig, Jr.\* The University of Texas at Austin Austin TX 78712  
+
+# Abstract  
+
+# 1. Introduction  
+
+Since the 1960s, substructuring, or component mode synthesis (CMS), has been used to model complex structures.  Substructuring involves dividing the structure into a number of substructures, or components (Fig. 1), obtaining reduced-order models of the components, and then assembling a reduced-order model of the entire structure. This paper defines and illustrates many of the terms that are found in the literature on substructuring, (fixed-interface modes, free-interface modes, constraint modes, residual fexibility, etc.), discusses a general procedure for coupling substructures, and compares two widely used methods. The focus of this paper is on the presentation of figures that illustrate the physical meaning of various component mode transformations.  
+
+![](4ad197aac2603d6a6f3df2a7b798d7ef0fb147f719dad4a067b5452e8cdd6d3d.jpg)  
+Figure 1: Typical Airplane Substructures.  
+
+When a large, complex structural system must be analyzed for its response to dynamic excitation, some form of substructure coupling method, or component mode synthesis (CMS) method, is usually employed. The term component modes is used to signify Ritz vectors, or assumed modes[1], that are used as basis vectors in describing the displacement of points within a substructure, or component. Component normal modes, or eigenvectors, are just one class of assumed modes. In the mid- $\cdot1960^{\circ}\mathrm{s}$ Hurty published several reports and papers on substructure coupling (e.g., [2, 3]). In collaboration with Hurty, Bamford created a CMS computer program that employed normal modes, rigid-body modes, constraint modes, and attachment modes[4]. A simplification of Hurty's method was presented by the author in 1968[5], and in the early 1970's MacNeal and Rubin introduced important alternatives to Hurty's CMS method[6, 7]. A number of CMS methods are described and compared in Refs. [8-10] and in at least three textbooks[11-13].  Although CMS methods have been developed for damped systems as well as for undamped systems, methods for damped systems are not discussed in the present paper.  
+
+Component mode synthesis involves three basic steps: division of a structure into components, definition of sets of component modes, and coupling of the component mode models to form a reduced-order system model. The primary uses of dynamic substructuring are: (1) to couple reduced-order models of moderately complex structures (e.g., airplane components, as in Fig. 1, or systems of automotive components), (2) in test-verification of finite element models of components, or (3) to implement computation of the dynamics of very large finite element models (e.g., multi-million-DOF models). This paper addresses primarily applications of the first type; Refs. [14--20] illustrate the relationship of substructure analysis to substructure testing, and Refs. [21, 22] are representative of the third application of substructuring.  
+
+In Section 2, the systematic procedures used to generate FE-based component modes are described, including discussions of inertia relief and residual fexibility. Section 3 is a review of a Lagrangemultiplier-based generalized substructure coupling procedure. This is followed, in Section 4, by a discussion of coupling analyses based on two widely used CMS methods, and, in Section 5, by conclusions.  
+
+# 2. Component Modes  
+
+The most general type of component, or substructure, is one that is connected to one or more adjacent components by redundant interfaces. Figure 2 illustrates a simple cantilever beam that is divided into three components; the middle one is a typical component with redundant interface (boundary) coordinates.  
+
+![](3afda16598240e86842e604d3ae7f3c497e9055c96b5327c4b7bcbb6ec181f01.jpg)  
+a. Components and coupled system.  
+
+![](973bf8ca2cdd7d1654e066b492c9c21debc986eb3758ee1e1a9a806b1682d9fe.jpg)  
+
+b. Typical component with redundant boundary.  
+
+As noted in Fig. 2, the coordinate sets $\tau,\mathcal{R},\mathcal{E}$ and $\pmb{{\cal B}}$ denote interior coordinates (i.e., not shared with an adjacent component), rigid-body coordinates, excess coordinates (i.e., redundant boundary coordinates), and boundary coordinates (i.e., shared with adjacent components). The numbers of coordinates in these sets are $N_{i},\,N_{r},\,N_{e},$ and $N_{b}$ , respectively, with $N_{b}=N_{r}+N_{e}$ and $N=N_{i}+N_{b}$  
+
+The equation of motion of $\mathbf{a}$ typical undamped component, labeled $^c$ may be written as  
+
+$$
+M^{c}\ddot{\boldsymbol{u}}^{c}+K^{c}\boldsymbol{u}^{c}=\boldsymbol{f}^{c}
+$$  
+
+where the superscript $c$ is the label of the particular component, and where $M^{c},\,K^{c}$ , and $\pmb{u}^{c}$ , are the component's mass matrix, stiffness matrix, and displacement vector, respectively. The force vector, $\pmb{f}^{c}$ includes both the externally applied forces and the forces on the component due to its connection to adjacent components at boundary degrees of freedom.  
+
+In component mode synthesis, the component's physical displacement coordinates $\pmb{\mathit{u}}$ are represented in terms of component generalized coordinates $\pmb{p}$ by the Rayleigh-Ritz coordinate transformation  
+
+$$
+\pmb{u}^{c}=\pmb{\Psi}^{c}\pmb{p}^{c}
+$$  
+
+where the component mode matrir $\Psi^{c}$ is a coordinate transformation matrix of preselected component (assumed) modes, including the following types: rigid-body modes, normal modes of free vibration (i.e., eigenvectors), constraint modes, and attachment modes. In collaboration with Hurty, Bamford defined all four of these types of modes in Ref. [4]; they are also defined in Refs. [9-11]. Other types of assumed modes (e.g., Krylov vectors [23]) may also be employed as component modes.  
+
+The coordinate transformation relating component physical coordinates ${\pmb u}^{c}$ to component generalized coordinates $\pmb{p}^{c}$ is given by Eq. (2). This equation, together with the equation of motion in generalized coordinates, forms the component modal model. From Eqs. (1) and (2), the component equation of motion in generalized coordinates is  
+
+$$
+\bar{M}^{\mathrm{c}}\,\ddot{\pmb{p}}^{\mathrm{c}}+\bar{K}^{\mathrm{c}}\,\pmb{p}^{\mathrm{c}}=\bar{\pmb{f}}^{c}
+$$  
+
+where  
+
+$$
+\bar{M}^{c}=\bar{\Psi}^{c T}M^{c}\Psi^{c},~\bar{K}^{c}=\bar{\Psi}^{c T}K^{c}\Psi^{c},~\bar{\pmb f}^{c}=\bar{\Psi}^{c T}\pmb f^{c}
+$$  
+
+The following partitioned forms of Eq. (1) will be useful in the derivation of component modes:  
+
+$$
+\begin{array}{r l r}&{}&{\left[\begin{array}{c c}{M_{i i}}&{M_{i b}}\\ {M_{b i}}&{M_{b b}}\end{array}\right]\left\{\begin{array}{c}{\ddot{u}_{i}}\\ {\ddot{u}_{b}}\end{array}\right\}+\left[\begin{array}{c c}{K_{i i}}&{K_{i b}}\\ {K_{b i}}&{K_{b b}}\end{array}\right]\left\{\begin{array}{c}{{u_{i}}}\\ {{u_{b}}}\end{array}\right\}}\\ &{}&{=\left\{\begin{array}{c}{{{f}_{i}}}\\ {{{f}_{b}}}\end{array}\right\}}\end{array}
+$$  
+
+$$
+\begin{array}{r l}&{\left[\begin{array}{c c c}{M_{\mathrm{ii}}}&{M_{\mathrm{ie}}}&{M_{\mathrm{ir}}}\\ {M_{\mathrm{ei}}}&{M_{\mathrm{ee}}}&{M_{\mathrm{er}}}\\ {M_{\mathrm{ri}}}&{M_{\mathrm{re}}}&{M_{\mathrm{rr}}}\end{array}\right]\left\{\begin{array}{c}{\dot{u}_{i}}\\ {\dot{u}_{e}}\\ {\dot{u}_{r}}\end{array}\right\}}\\ &{+\left[\begin{array}{c c c}{K_{\mathrm{ii}}}&{K_{\mathrm{ie}}}&{K_{\mathrm{ir}}}\\ {K_{\mathrm{ei}}}&{K_{\mathrm{ee}}}&{K_{\mathrm{er}}}\\ {K_{\mathrm{ri}}}&{K_{\mathrm{re}}}&{K_{\mathrm{rr}}}\end{array}\right]\left\{\begin{array}{c}{u_{i}}\\ {u_{e}}\\ {u_{r}}\end{array}\right\}}\\ &{=\left\{\begin{array}{c}{f_{i}}\\ {f_{e}}\\ {f_{r}}\end{array}\right\}}\end{array}
+$$  
+
+The superscript $c$ ,which was used above to designate a component, will be omitted from component matrices and vectors in the remainder of this section.  
+
+# 2.1. Normal Modes  
+
+Component normal modes are eigenvectors, and may be classified according to the interface boundary conditions specified for the component - fixed-interface normal modes, free-interface normal modes, hybrid-interface normal modes, or loadedinterface  normal  modes[11].  Component  ficedinterface normal modes are obtained by restraining all boundary DOFs and solving the following eigenproblem:  
+
+$$
+\left[K_{i i}-\omega_{j}^{2}M_{i i}\right]\left\{\phi_{i}\right\}_{j}=0,\;j=1,2,...,N_{i}
+$$  
+
+The complete set of $N_{i}$ fixed-interface normal modes is labeled $\Phi_{n}$ and assembled according to the partitioning of Eq. (5) as columns of the modal matrix  
+
+$$
+\boldsymbol{\Phi}_{n}\equiv\left[\begin{array}{l}{\boldsymbol{\Phi}_{i n}}\\ {\boldsymbol{0}_{b n}}\end{array}\right]
+$$  
+
+When normalized with respect to the mass matrix $M_{i i}$ , the fixed-interfacc modcs satisfy  
+
+$$
+\Phi_{i n}^{T}M_{i i}\Phi_{i n}=I_{n n},\ \Phi_{i n}^{T}K_{i i}\Phi_{i n}=\Lambda_{n n}=\mathrm{diag}(\omega_{j}^{2})
+$$  
+
+Figure 3 shows the $N_{i}$ (four) fixed-interface normal modes for the 8DOF free-free beam component in Fig. 2b.  
+
+![](ddfd148b5e4c3e6a85e278be7529012c90ffed8461442fd5f6920dd8fcb54646.jpg)  
+Figure 3: Component Fixed-Interface Modes  
+
+A second type of component normal modes used in CMS is free-interface normal modes. These normal modes are defined by the equation  
+
+$$
+\left[K-\omega_{j}^{2}M\right]\{\phi\}_{j}=0,\;j=1,2,...,(N_{f}=N-N_{r}),
+$$  
+
+The assembled set of $N_{f}$ flexible (i.e., non rigidbody) free-interface normal modes is  
+
+$$
+\boldsymbol{\Phi}_{n}\equiv\left[\begin{array}{l}{\boldsymbol{\Phi}_{i n}}\\ {\boldsymbol{\Phi}_{b n}}\end{array}\right]
+$$  
+
+Figure 4 shows the first four of the six free-free fex modes of the 8DOF component.  
+
+![](ca543703826cffbe53c17075f16e1e82d1836850b2d8a5b7a4a3e55e9281e91d.jpg)  
+Figure 4: Flex Modes for the Free-Free Component.  
+
+A third important type of component normal modes is loaded-interface normal modes. This includes lumped-mass loaded-interface normal modes, commonly referred to as mass-additive normal modes[16, 19]. Benfield and Hruda described CMS methods based on “consistent" mass- and stiffnessadditive normal modes[24], but these methods require reduced-order models of all adjacent substructures, so they are generally of limited practical value.  
+
+# 2.2. Constraint Modes; Rigid-body Modes  
+
+A constraint mode is defined as the static deformation of a structure when a unit displacement is applied to one coordinate of a specified set of “constraint" coordinates, $c$ , while the remaining coordinates of that set are restrained, and the remaining degrees of freedom of the structure are force-free. The set of interface constraint modes based on unit displacement of the boundary coordinates $\pmb{u}_{b}$ is a very useful CMS set, because of the ease of enforcing inter-component compatibility when these constraint modes are employed, as will be explained in Section 4.1. This set, with ${\mathcal{C}}={\boldsymbol{B}}$ , is given by  
+
+$$
+\left[\begin{array}{c c}{{K_{i i}}}&{{K_{i b}}}\\ {{K_{b i}}}&{{K_{b b}}}\end{array}\right]\left[\begin{array}{c}{{\Psi_{i b}}}\\ {{I_{b b}}}\end{array}\right]=\left[\begin{array}{c}{{0_{i b}}}\\ {{R_{b b}}}\end{array}\right]
+$$  
+
+That is, the constraint-mode matrix $\Psi_{c}$ is given by  
+
+$$
+\boldsymbol{\Psi}_{c}\equiv\left[\begin{array}{c}{\boldsymbol{\Psi}_{i b}}\\ {\boldsymbol{I}_{b b}}\end{array}\right]=\left[\begin{array}{c}{-\boldsymbol{K}_{i i}^{-1}\boldsymbol{K}_{i b}}\\ {\boldsymbol{I}_{b b}}\end{array}\right]
+$$  
+
+From Eqs. (8) and (12) it can easily be shown that these constraint modes are stiffness-orthogonal to all of the fixed-interface normal modes, that is,  
+
+$$
+\Phi_{n}^{T}K\bar{\Psi}_{c}=0
+$$  
+
+Figure 5 shows the $N_{b}$ (four) constraint modes for the 8DOF free-free beam component in Fig. 2b.  
+
+![](3be0229b946c4efb71ed5d9586cba4305d6b05d7395a4cc545b3cae3c56c9f55.jpg)  
+Figure 5: Component Constraint Modes.  
+
+Although they are often considered to be normal modes, rigid-body modes are actually a special case of constraint modes. They can be defined relative to anysetof $N_{r}$ coordinates that is just sufficient to restrain rigid-body motion of the component. For purposes of substructure coupling, rigid-body modes will be defined relative to a set $\mathcal{R}$ of boundary coordinates. Then,  
+
+$$
+\left[\begin{array}{l l l}{K_{i i}}&{K_{i e}}&{K_{i r}}\\ {K_{e i}}&{K_{e e}}&{K_{e r}}\\ {K_{r i}}&{K_{r e}}&{K_{r r}}\end{array}\right]\left[\begin{array}{l}{\Psi_{i r}}\\ {\Psi_{e r}}\\ {I_{r r}}\end{array}\right]=\left[\begin{array}{l}{0_{i r}}\\ {0_{e r}}\\ {0_{r r}}\end{array}\right]
+$$  
+
+so the set of rigid-body modes is obtained by solving the top two row partitions of Eq. (15), giving  
+
+$$
+\boldsymbol{\Psi}_{r}\equiv\left[\begin{array}{l}{\Psi_{i r}}\\ {\Psi_{e r}}\\ {I_{r r}}\end{array}\right]=\left[\begin{array}{l l}{-\left[\begin{array}{l l}{G_{i i}}&{G_{i e}}\\ {G_{e i}}&{G_{e e}}\end{array}\right]\left[\begin{array}{l}{K_{i r}}\\ {K_{e r}}\end{array}\right]}\\ {I_{r r}}\end{array}\right]_{\ldots}
+$$  
+
+where  
+
+$$
+\bar{G}_{c}\equiv\left[\begin{array}{l l}{G_{i i}}&{G_{i e}}\\ {G_{e i}}&{G_{e e}}\end{array}\right]=\left[\begin{array}{l l}{K_{i i}}&{K_{i e}}\\ {K_{e i}}&{K_{e e}}\end{array}\right]^{-1}
+$$  
+
+is the cantilever fleribility matriz for the component restrained at the $\mathcal{R}$ coordinates.  
+
+Redundant-interface constraint modes can then be defined for unit displacements at the redundant  
+
+(excess) boundary coordinate set $\mathcal{E}$ , and with the $\mathcal{R}$ coordinates fixed, by the equation  
+
+$$
+\left[\begin{array}{c c c}{K_{i i}}&{K_{i e}}&{K_{i r}}\\ {K_{e i}}&{K_{e e}}&{K_{e r}}\\ {K_{r i}}&{K_{r e}}&{K_{r r}}\end{array}\right]\quad\left[\begin{array}{c}{\Psi_{i e}}\\ {I_{e e}}\\ {0_{r e}}\end{array}\right]=\left[\begin{array}{c}{0_{i e}}\\ {R_{e e}}\\ {R_{r e}}\end{array}\right]
+$$  
+
+Then,  
+
+$$
+\underbrace{\Psi_{e}}_{N\times N_{e}}\equiv\left[\begin{array}{c}{\Psi_{i e}}\\ {I_{e e}}\\ {0_{r e}}\end{array}\right]=\left[\begin{array}{c}{-K_{i i}^{-1}K_{i e}}\\ {I_{e e}}\\ {0_{r e}}\end{array}\right]
+$$  
+
+Either the set of interface constraint modes $\Psi_{c}$ defined by Eq. (13), or the combined set $\left[\Psi_{r}\ \Psi_{e}\right]$ defined by Eqs. (16) and (19), spans the static response of the substructure to interface loading and allows for arbitrary interface displacements $\pmb{u}_{b}$ . Along with the interface displacement, there is accompanying displacement of the interior of the substructure, as determined by Eqs. (13), (16), and (19). Additional interior fexibility can be incorporated by including fixed-interface normal modes, fixed-interface Krylov vectors, or other fixed-interface assumed modes in the component mode matrix $\Psi[3,\,5,\,23]$  
+
+# 2.3. Attachment Modes  
+
+An attachment mode is defined as the component displacement vector due to a single unit force applied at one of the coordinates of a given set $\pmb{A}$ Consequently, attachment modes are just columns of the associated fexibility matrix. Attachment modes were defined by Bamford[4], and they get their name from their usefulness in representating the deformation of a structure to loading (e.g., an external force, an attached mass, or an attached fexible component) at the point where the attachment mode's unit force is applied. In this paper we are interested in defining attachment modes to represent the response of a component to forces at its interface with adjoining components. One diffculty encountered in using attachment modes is that many components have one to six rigid-body degrees of freedom, making it impossible to apply directly to the unrestrained component the necessary unit forces in order to compute the resulting attachment mode shapes. However, one option in this case is to select a set $\mathcal{R}$ of boundary rigid-body degrees of freedom, (mathematically) restrain the component at these DOFs, and then form cantilever attachment modes by applying unit loads at the redundant boundary coordinates, that is, for $\boldsymbol{A}=\boldsymbol{\mathcal{E}}$ .Then,  
+
+$$
+\left[\begin{array}{l l l}{K_{i i}}&{K_{i e}}&{K_{i r}}\\ {K_{e i}}&{K_{e e}}&{K_{e r}}\\ {K_{r i}}&{K_{r e}}&{K_{r r}}\end{array}\right]\left[\begin{array}{l}{\Psi_{i e}}\\ {\Psi_{e e}}\\ {0_{r e}}\end{array}\right]=\left[\begin{array}{l}{0_{i e}}\\ {I_{e e}}\\ {R_{r e}}\end{array}\right]
+$$  
+
+It can be seen that these attachment modes are just an expanded form of the columns of the righthand partition of the Hexibility matrix $\bar{G}_{c}$ of Eq. (17) with $A=\mathcal{E}$ . That is,  
+
+$$
+\boldsymbol{\Psi}_{\boldsymbol{s}\,\boldsymbol{N}_{c}}\equiv\left[\begin{array}{l}{\Psi_{i e}}\\ {\Psi_{e e}}\\ {0_{r e}}\end{array}\right]=\left[\begin{array}{l}{G_{i e}}\\ {G_{e e}}\\ {0_{r e}}\end{array}\right]
+$$  
+
+Two important topics that arise when freeinterface normal modes are to be employed to represent the fexible behavior of unrestrained components are inertia relief and residual fezibility, both of which were discussed by MacNeal[6] and Rubin[7]. Sections 2.4 and 2.5 treat these two topics, and the related forms of attachment modes are defined.  
+
+# 2.4. Inertia-Relief Attachment Modes  
+
+When a component has rigid-body freedom, it is appropriate to employ inertia-relief attachment modes[6, 7, 11]. The term inertia relief refers to the process of applying to the component an equilibrated load system $\pmb{f}_{f}$ , which consists of the original force vector $\pmb{f}$ equilibrated by the rigid-body d'Alembert force vector $M\ddot{\boldsymbol{u}}_{r}$ , where $\pmb{u}_{r}$ is the rigid-body motion due to $\pmb{f}$ .Starting with Eq. (1), let the displacement vector be separated into rigid-body displacement and fexible-body displacement, that is, let  
+
+$$
+\pmb{u}=\pmb{u}_{r}+\pmb{u}_{f}=\Psi_{r}\pmb{p}_{r}+\Phi_{f}\pmb{p}_{f}
+$$  
+
+whcrc all of thc $N_{f}$ fexiblc-body modes are included in $\Phi_{f}$ . Then, the equations  
+
+$$
+\Psi_{r}^{T}M\Phi_{f}=0\ ,\ \ \bar{M}_{r r}=\Psi_{r}^{T}M\Psi_{r}
+$$  
+
+are the appropriate orthogonality equation and the definition of the rigid-body modal mass matrix, respectively. (It is not assumed that the rigid-body modes are orthonormalized.） Since $K\Psi_{r}~\,=\,\,0$ Eqs. (1) and (22) can be combined to give  
+
+$$
+M\Phi_{f}\ddot{p}_{f}+K\Phi_{f}p_{f}=f-M\Psi_{r}\ddot{p}_{r}
+$$  
+
+When this equation is premultiplied by $\Psi_{r}^{T}$ and orthogonality is invoked, the result is the selfequilibrated force  
+
+$$
+\pmb{f}_{f}=\pmb{f}-M\ddot{\pmb{u}}_{r}=P_{r}\pmb{f}
+$$  
+
+where $P_{r}$ is the inertia-relief projection matriz, de fined by  
+
+$$
+P_{r}=I-M\Psi_{r}\bar{M}_{r r}^{-1}\Psi_{r}^{T}
+$$  
+
+When any force vector is premultiplied by this inertia-relief projection matrix, the corresponding  
+
+force system is self-equilibrated. Also, from Eq. (26) it can easily be verified that $P_{r}^{T}$ is mass-orthogonal to the rigid-body modes, that is,  
+
+$$
+\Psi_{r}{}^{T}M P_{r}^{T}=0
+$$  
+
+Inertia-relief attachment modes are staticdeformation shapes defined by applying unit forces at the all interface coordinates $\quad A=B!$ , that is, by applying the force  
+
+$$
+F_{b}=\left[\begin{array}{c}{{0_{i b}}}\\ {{I_{b b}}}\end{array}\right]
+$$  
+
+pre-multiplied by the inertia-relief projection matrix, $P_{r}$ . Since the unit-force column vectors in $F_{b}$ are self-equilibrated by the inertia-relief projection matrix, no reaction forces are required, such as there are in Eq. (20). Deformation of the component due to this equilibrated force system is given by  
+
+$$
+\hat{\Psi}_{b}=G_{c}P_{r}F_{b}
+$$  
+
+where $G_{c}$ is the constrained feribility matriz, a spe cial expanded (singular) form of the cantilever fexibility matrix $\dot{G}_{c}$ in Eq. (17), given by  
+
+$$
+G_{c}={\left[\begin{array}{l l l}{G_{i i}}&{G_{i e}}&{0_{i r}}\\ {G_{e i}}&{G_{e e}}&{0_{e r}}\\ {0_{r i}}&{0_{r e}}&{0_{r r}}\end{array}\right]}
+$$  
+
+The attachment-mode set defined by Eq. (29) is made orthogonal to the rigid-body modes, and the resulting inertia-relief attachment modes are given by[11]  
+
+$$
+\Psi_{b}\equiv\left[\begin{array}{c}{{\Psi_{i b}}}\\ {{\Psi_{b b}}}\end{array}\right]=G_{f}\left[\begin{array}{c}{{0_{i b}}}\\ {{I_{b b}}}\end{array}\right]
+$$  
+
+where  
+
+$$
+G_{f}=P_{r}^{T}G_{c}P_{r}
+$$  
+
+is the elastic flezibility matrir in inertia-relief format. In Eq. (31), the $\boldsymbol{I_{b b}}$ matrix in $F_{b}$ picks out the columns of the fexibility matrix $G_{f}$ that correspond to unit forces applied at the boundary. From Eqs. (27) and (32), it can be shown that the columns of $G_{f}$ are mass-orthogonal to the rigid-body modes $\Psi_{r}$ . Therefore, $G_{f}$ spans the same subspace as do the free-interface fex modes of Eq. (11).  
+
+The top two plots in Fig. 6 are the shapes that correspond to the two columns of the elastic flexibility matrix $G_{f}$ for unit forces at the transverse DOFs at the two ends of the 8DOF free-free beam in Fig. 2. It is clear that these two fexibility shapes (and the remaining six as well) are dominated by the contribution of the fundamental free-free fex mode (Fig. 4). This can be seen clearly by the bottom two figures, which represent symmetric and antisymmetric loading by unit forces at the two ends of the component.  
+
+![](e2bae7628e89c7f6bfb5aea5373825281f97d837709c0f407c22637999b7a52a.jpg)  
+Figure 6: Elastic Flexibility Shapes  
+
+# 2.5. Flexibility Matrices and Residual Flexibility Attachment Modes  
+
+The complete set of component normal modes $\Phi_{n}$ and the corresponding set of eigenvalues $\Lambda_{n n}$ are identified by the subscript $n$ , whether these are the $N_{i}$ fixed-interface modes, the $N_{f}$ free-free flexible (fex) modes, or some other form of component normal modes.  
+
+Let the (diagonal) modal mass matrix and modal stiffness matrix for modes $\Phi_{n}$ be  
+
+$$
+{\bar{M}}=\Phi_{n}^{T}M\Phi_{n}\ ,\ \ {\bar{K}}=\Phi_{n}^{T}K\Phi_{n}
+$$  
+
+respectively. Then, the elastic feribility matrir, $G$ of the component can be expressed in the following mode-superposition format:  
+
+$$
+G=\Phi_{n}\bar{K}^{-1}\Phi_{n}^{T}=\sum_{j}\phi_{j}\left(\frac{1}{\bar{K}_{j}}\right)\phi_{j}^{T}
+$$  
+
+Note that each column of the $j$ th mode's contribution to the elastic fexibility matrix has the shape of mode $\phi_{j}$ .Although the elastic fexibility matrices $G$ of Eq. (34) and matrix $G_{f}$ of Eq. (32) and illustrated in Fig. 6 are formed in different ways, they are numerically the same. In this section we will be concerned with components that have rigid-body freedom, in which case $G$ is singular, with rank $N_{f}$ Regardless of whether the elastic fexibility matrix  
+
+$G$ is singular or not, from Eqs. (33b) and (34) it can be shown that  
+
+$$
+G^{T}K G=G
+$$  
+
+Since model reduction is one of the major objectives in CMS, the normal mode set is usually reduced to a smaller set of kept normal modes, denoted by $\Phi_{k}$ ，where $\Phi_{n}~\equiv~\left[\Phi_{k}~\Phi_{d}\right]$ .\*The deleted normal modes, $\Phi_{d}$ , are generally all of the modes above some specified cutof frequency. The portion of the fexibility matrix contributed by modes $\Phi_{d}$ is called the residual feribility matrir. It is given by  
+
+$$
+G_{d}=\Phi_{d}\bar{K}_{d d}^{-1}\Phi_{d}^{T}=G-\Phi_{k}\bar{K}_{k k}^{-1}\Phi_{k}^{T}
+$$  
+
+where $G$ is the total fexibility matrix. Since it is not usually feasible to compute or measure the $\Phi_{d}$ modes, Eq. (36) is useful only because Eq. (32) exists as an alternative to Eq. (34) for determining the elastic fexibility matrix $G$  
+
+The matrix $G_{d}$ will always be a singular matrix because of the modes deleted in Eq. (36). Also, because of the mass- and stiffness-orthogonality of the kept modes to the deleted modes,  
+
+$$
+\Phi_{k}^{T}M G_{d}=0,\ \mathrm{and}\ \ \Phi_{k}^{T}K G_{d}=0
+$$  
+
+and because of the orthogonality between all rigidbody modes and all fexible-body modes,  
+
+$$
+\Psi_{r}^{T}M G_{d}=0,\ \mathrm{and}\ \ \Psi_{r}^{T}K G_{d}=0
+$$  
+
+Residual-fleribility attachment modes maybe de fined for forces applied at the interface coordinates, that is, for $\boldsymbol{A}=\boldsymbol{B}$ , by the following equation:  
+
+$$
+\begin{array}{r c l}{{\Psi_{d}}}&{{\equiv}}&{{\left[\begin{array}{c}{{\Psi_{i b}}}\\ {{\Psi_{b b}}}\end{array}\right]=G_{d}F_{b}}}\\ {{}}&{{=}}&{{[G_{f}-\Phi_{k}\bar{K}_{k k}^{-1}\Phi_{k}^{T}]\left[\begin{array}{c}{{0_{i b}}}\\ {{I_{b b}}}\end{array}\right]}}\end{array}
+$$  
+
+Figure 7 shows the attachment mode shape for the component with a unit force at DOF7; the top figure includes all six flex modes, the middle figure is the contribution of two “kept" modes, and the bottom figure is the corresponding residual fexibility attachment mode shape. It is clear that the order of magnitude of the residual fexibility is smaller than that of the fexibility of the kept modes. Figure 8 shows the residual-Hexibility attachmentmode shapes $({\bf k}{=}2)$ for the component with unit forces at DOFs 7 and 5 (left and right ends). The top two figure are the attachment-mode shapes for the individual unit forces; the middle figure is the shape produces by symmetric loading by two unit forces, and the bottom figure is the corresponding residual-fexibility attachment-mode shape for antisymmetric loading. It can easily be seen that these residual-fexibility shapes are free of the first two (kept) normal-mode contributions.  
+
+![](88ff35243039fc2e46d94f50e1031720e313dca7bbafb0eba9b0054bc93ed7b5.jpg)  
+Figure 7: An Illustration of Residual Flexibility.  
+
+![](7792deafc9126c656d99fceec86135d06a19a7955648c01f9bb38c89e622fa62.jpg)  
+Figure 8: More Residual-Flexibility Shapes.  
+
+Incorporation of $\Psi_{d}$ into the component mode set ensures complete representation of static defection of the component due to forces applied at interface DOFs. In this sense, it is closely related to the mode-acceleration method for incorporating static completeness in dynamic-response computations[1, 7, 11]. Hintz has given an extensive discussion of the need for statically complete component mode sets in Ref. [25]. We will return to the topic of residual flexibility later in Section 4.2.  
+
+# 3. A Generalized Component Coupling Procedure for Undamped Structures  
+
+In this section a generalized substructure coupling procedure that employs Lagrange multipliers to enforce  inter-component  displacement compatibility equations (and other constraint equations, if applicable) is presented. Let the system be composed of two components, labeled $_\alpha$ and $\beta$ , that have a common (generally redundant) interface. The physical displacements at the interface are constrained by the displacement compatibility equation  
+
+$$
+\pmb{u}_{b}^{\alpha}=\pmb{u}_{b}^{\beta}
+$$  
+
+and the mutually reactive interface forces (i.e., not including external forces applied at the interface) are related by  
+
+$$
+\hat{\pmb f}_{b}^{\alpha}+\hat{\pmb f}_{b}^{\beta}=\mathbf{0}
+$$  
+
+Constraint equations, such as Eq. (40) and any other constraint equations that are to be imposed (say $N_{C}$ equations in all), can be written in terms of the generalized coordinates $\pmb{p}$ and combined to form a matrix constraint equation of the form  
+
+$$
+C{\pmb p}={\bf0}
+$$  
+
+For example, Eqs. (2) and (40) can be combined to give the constraint equation  
+
+$$
+\left[\Psi_{b}^{\alpha}\;-\Psi_{b}^{\beta}\right]\left\{\begin{array}{l}{{p^{\alpha}}}\\ {{p^{\beta}}}\end{array}\right\}={\bf0}
+$$  
+
+The synthesis of the system equation of motion is based on Lagrange's equation of motion with undetermined multipliers[11, 26]. The Lagrangian for the system of two coupled substructures can be written  
+
+$$
+\mathcal{L}=\mathcal{T}-\mathcal{V}+\lambda^{T}C p
+$$  
+
+where $\tau$ is the system kinetic energy and $\nu$ is the system potential energy, given by  
+
+$$
+\begin{array}{l}{{{\cal T}=\frac{1}{2}\dot{\pmb{p}}^{\alpha T}\bar{M}^{\alpha}\dot{\pmb{p}}^{\alpha}+\frac{1}{2}\dot{\pmb{p}}^{\beta T}\bar{M}^{\beta}\dot{\pmb{p}}^{\beta}=\frac{1}{2}\dot{\pmb{p}}^{T}\bar{M}\dot{\pmb{p}}}}\\ {{\mathcal{V}=\frac{1}{2}\pmb{p}^{\alpha T}\bar{K}^{\alpha}\pmb{p}^{\alpha}+\frac{1}{2}\pmb{p}^{\beta T}\bar{K}^{\beta}\pmb{p}^{\beta}=\frac{1}{2}\pmb{p}^{T}\bar{K}\pmb{p}}}\end{array}
+$$  
+
+where  
+
+$$
+\begin{array}{r l}&{\bar{M}\equiv\left[\begin{array}{c c}{\bar{M}^{\alpha}}&{0}\\ {0}&{\bar{M}^{\beta}}\end{array}\right],~\bar{K}\equiv\left[\begin{array}{c c}{\bar{K}^{\alpha}}&{0}\\ {0}&{\bar{K}^{\beta}}\end{array}\right],}\\ &{p\equiv\left\{\begin{array}{c}{p^{\alpha}}\\ {p^{\beta}}\end{array}\right\}}\end{array}
+$$  
+
+Corresponding to Eq. (46),  
+
+$$
+\bar{\pmb{f}}\equiv\left\{\begin{array}{l}{\bar{\pmb{f}}^{\alpha}}\\ {\bar{\pmb{f}}^{\beta}}\end{array}\right\}=\left\{\begin{array}{l}{\Psi^{\alpha T}\pmb{f}^{\alpha}}\\ {\Psi^{\beta T}\pmb{f}^{\beta}}\end{array}\right\}
+$$  
+
+The system equations of motion can now be obtained by applying Lagrange's equation in the form  
+
+$$
+\frac{d}{d t}\left(\frac{\partial\mathcal{L}}{\partial\dot{p}_{j}}\right)-\frac{\partial\mathcal{L}}{\partial p_{j}}=\bar{f}_{j}\;\;,\;\;\;j=1,2,...,N_{\alpha}\,+N_{\beta}
+$$  
+
+where $p_{j}$ refers to the $j$ th element of the merged displacement vector $\pmb{p}$ , and ${\bar{f}}_{j}$ refers to the corre sponding (externally applied) force. As required by Eq. (41), the mutually reactive interface forces cancel out and do not appear on the right-hand side of Eq. (48). In matrix form, the $\left(N_{\alpha}+N_{\beta}\right)$ equations of Eq. (48) can be written as  
+
+$$
+\bar{M}\ddot{\pmb{p}}+\bar{K}\pmb{p}=\bar{\pmb{f}}+C^{T}\pmb{\lambda}
+$$  
+
+Since, due to the constraint equation, Eq. (42), thecoordinates $\pmb{p}$ are not linearly independent, practically all substructure coupling methods solve the coupled set of equations, Eqs. (42) and (49), by introducing a linear transformation of the form  
+
+$$
+\pmb{p}=S\pmb{q}
+$$  
+
+where $\pmb q$ is the vector of independent system generalized coordinates.  
+
+Let $\pmb{p}$ be rearranged, if necessary, and partitioned into $N_{C}$ dependent coordinates $\pmb{p}_{D}$ , and $(N_{\alpha}\!+\!N_{\beta}-$ $N_{C}\,.$ ) linearly independent coordinates $\pmb{p}_{I}$ , and let Eq. (42) be partitioned accordingly, giving  
+
+$$
+\left\{C_{\scriptscriptstyle D\scriptscriptstyle D}C_{\scriptscriptstyle D\scriptscriptstyle I}\right\}\left\{\begin{array}{c}{{p_{\scriptscriptstyle I\!}}}\\ {{p_{\scriptscriptstyle I\!}}}\end{array}\right\}={\bf0}
+$$  
+
+where $C_{D D}$ is a nonsingular square matrix. Then, the equation  
+
+$$
+p\equiv\left\{\begin{array}{c}{{p_{n}}}\\ {{p_{I}}}\end{array}\right\}=\left[\begin{array}{c}{{-C_{\scriptscriptstyle D D}^{-1}C_{\scriptscriptstyle D I}}}\\ {{I_{\scriptscriptstyle I I}}}\end{array}\right]p_{I}\equiv S q
+$$  
+
+defines both $S$ and $\pmb q$ . Then, the vector of independent system generalized coordinates is $\pmb q\equiv\pmb{p}_{I}$ ,and the coupling transformation matrix $S$ is given by  
+
+$$
+S=\left[\begin{array}{c}{{-C_{{}_{D D}}^{-1}C_{{}_{D I}}}}\\ {{{}_{I_{I I}}}}\end{array}\right]
+$$  
+
+Substitution of Eq. (50) into Eq. (49) and premultiplication of the resulting equation by $S^{T}$ gives  
+
+$$
+\begin{array}{r}{M_{q}\ddot{\pmb q}+K_{q}\pmb q=\pmb f_{q}+\pmb S^{T}\pmb C^{T}\pmb\lambda}\end{array}
+$$  
+
+where  
+
+$$
+M_{q}=S^{T}\bar{M}S,\ K_{q}=S^{T}\bar{K}S,\ f_{q}=S^{T}\bar{f}
+$$  
+
+From Eqs. (51) and (52), it is seen that $C S\,=\,0$ Therefore, the system equation of motion, Eq. (54), becomes simply  
+
+$$
+M_{q}{\ddot{q}}+K_{q}q=f_{q}
+$$  
+
+Although Eq.(55) defines $M_{q},~K_{q}$ ,and $\pmb{f}_{q}$ in terms of matrix operations, the system matrices and force vector can usually be assembled from the substructure matrices by the “direct stiffness" assembly procedure, as is illustrated in Section 4.1.  
+
+Equations (42) through (56) describe a single level of substructuring; however, essentially the same procedure can be employed when a structure is partitioned into several levels of substructures (e.g.. Ref.[21]).  
+
+# 4. Component Mode Synthesis Methods  
+
+Most applications of component mode synthesis employ one of two approaches, called constraint-mode methods and attachment-mode methods. The former employ constraint modes and fixed-interface normal modes, as represented by Hurty's method[3] and the Craig-Bampton variant of Hurty's method[5]. The latter employ attachment modes and freeinterface normal modes, as represented by MacNeal's method[6] and Rubin's method[7]. It is possible to cite here only a small number of the significant papers dealing with the use of component modes in structural dynamics.  
+
+# 4.1 Constraint-Mode Methods  
+
+Although there had been previous applications of component modes, Hurty's 1965 paper[3] provided the first comprehensive development of a finite element oriented CMS method based on constraint modes and fixed-interface modes. Craig and Bampton[5] simplified Hurty's method by treating all interface degrees of freedom together, rather than requiring the interface degrees of freedom to be separated into rigid-body freedoms and redundant interface freedoms. The displacement transformation for this method employs a combination of fixed-interface normal modes (Eq. (7) and Fig. 3) and interface constraint modes (Eq. (13) and Fig. 5), and takes the form  
+
+$$
+\begin{array}{r}{\boldsymbol{u}^{c}\equiv\left\{\begin{array}{c}{\boldsymbol{u}_{i}}\\ {\boldsymbol{u}_{b}}\end{array}\right\}^{c}=\left[\begin{array}{c c}{\boldsymbol{\Phi}_{i k}}&{\boldsymbol{\Psi}_{i b}}\\ {\boldsymbol{0}}&{I_{b b}}\end{array}\right]^{c}\left\{\begin{array}{c}{\boldsymbol{p}_{k}}\\ {\boldsymbol{p}_{b}}\end{array}\right\}^{c}}\end{array}
+$$  
+
+where the "Craig-Bampton transformation matrix" is  
+
+$$
+\Psi_{C B}^{c}=\left[\begin{array}{c c}{{\Phi_{i k}}}&{{\Psi_{i b}}}\\ {{0}}&{{I_{b b}}}\end{array}\right]^{c}
+$$  
+
+where $\Phi_{i k}$ is the interior partition of the fixedinterface modal matrix, and $\Psi_{i b}$ is the interior partition of the constraint-mode matrix.  
+
+With component fixed-interface normal modes normalized according to Eq. (9), the reduced mass and stiffness matrices, Eqs. (4), have the special forms  
+
+$$
+\begin{array}{r}{\bar{M}_{C B}^{c}=\left[\begin{array}{c c}{I_{k k}}&{\bar{M}_{k b}}\\ {\bar{M}_{b k}}&{\bar{M}_{b b}}\end{array}\right]^{c}\,,\;\;\bar{K}_{C B}^{c}=\left[\begin{array}{c c}{\Lambda_{k k}}&{0_{k b}}\\ {0_{b k}}&{\bar{K}_{b b}}\end{array}\right]^{c}}\end{array}
+$$  
+
+The zeros in the $k b$ and $b k$ partitions of $\bar{K}_{C B}^{c}$ arethe result of orthogonality equation Eq. (14).  
+
+Equation (57) implies that $\pmb{p}_{b}^{c}=\pmb{u}_{b}^{c}$ .Therefore, in terms of component generalized coordinates, the interface compatibility equation, Eq. (40), becomes  
+
+$$
+\pmb{p}_{b}^{\alpha}=\pmb{p}_{b}^{\beta}=\pmb{q}_{b}=\pmb{u}_{b}
+$$  
+
+Then, Eq. (50) takes the form  
+
+$$
+\left\{\begin{array}{c}{p_{k}^{\alpha}}\\ {p_{b}^{\alpha}}\\ {p_{k}^{\beta}}\\ {p_{b}^{\beta}}\end{array}\right\}=\left[\begin{array}{c c c}{I}&{0}&{0}\\ {0}&{0}&{I}\\ {0}&{I}&{0}\\ {0}&{0}&{I}\end{array}\right]\left\{\begin{array}{c}{q_{k}^{\alpha}}\\ {q_{k}^{\beta}}\\ {u_{b}}\end{array}\right\}
+$$  
+
+and the component coupling matrix $S$ is justthe "direct-stiffness assembly” matrix. For the twocomponent system, component mass and stiffness matrices (Eq. (59)) are assembled to form the following system reduced-order mass and stiffness matrices, respectively:  
+
+$$
+\begin{array}{r l}&{M_{C B}=\left[\begin{array}{l l l}{I_{k_{a}k_{o}}^{\alpha}}&{0_{k_{o}k_{i}\beta}^{\alpha}}&{\hat{M}_{k_{o}k\beta}^{\alpha}}\\ {0_{k_{i}k_{o}}^{\beta}}&{I_{k_{i}k_{i}\beta}^{\beta}}&{\hat{M}_{k_{i}\beta}^{\beta}}\\ {\bar{M}_{b k_{o}}^{\alpha}}&{\bar{M}_{b k_{o}}^{\beta}}&{\bar{M}_{b k}^{\alpha}+\bar{M}_{b k}^{\beta}}\end{array}\right]}\\ &{K_{C B}=\left[\begin{array}{l l l}{\Lambda_{k_{c}k_{o}}^{\alpha}}&{0_{k_{c}k_{i}}^{\alpha}}&{0_{k_{c}b}^{\alpha}}\\ {0_{k_{i}k_{o}}^{\beta}}&{\Lambda_{k_{i}k_{\beta}}^{\beta}}&{0_{k_{i}b}^{\beta}}\\ {0_{k_{k_{o}}}^{\alpha}}&{0_{k_{k_{\beta}}}^{\beta}}&{\bar{K}_{b k}^{\alpha}+\bar{K}_{b k}^{\beta}}\end{array}\right]}\end{array}
+$$  
+
+Thus, component models based on the use of fixed-interface modes plus interface constraint modes are essentially superelements - all physical boundary coordinates are retained as independent generalized coordinates, greatly facilitating component coupling. Because of the simple, straightforward procedures for formulating the component modes employed by this method, because of the straightforward way in which components are coupled to form the component mode system model and the sparsity patterns of the resulting system matrices, and because this method also produces highly accurate models with relatively few component modes[8], this method has been widely used and is available in a number of commercial finite element codes (e.g., MSC/NASTRAN[27]).  
+
+# 4.2 Attachment-Mode Methods  
+
+If a reduced set of component normal modes is used without including a complete set of either interface constraint modes or interface attachment modes, the component mode set is not statically complete, as indicated in Section 2.5. This is true of the “classical" CMS method of using only a set of free-interface normal modes[8]. However, methods that employ free-interface normal modes together with attachment modes (including residualfexibility attachment modes and/or inertia-relief attachment modes) are widely used, especially MacNeal's method and Rubin's methodl6, 7, 26], and especially in context of experimental verification of finite element models[14-20]. Martinez, et al.[14, 28], simplified the formulation of this class of CMS methods, casting the component mode matrix $\Psi$ inthe same format as that of the constraint-mode method of Craig and Bampton.  The variant of Rubin's method proposed by Martinez, et al., is described below.  
+
+The initial form of the displacement transformation for this method employs a combination of rigidbody modes (from Eq. (16)), kept free-free normal modes (from Eq. (11)), and residual-fexibility attachment modes (from Eq. (39)). To simplify the following derivation, the rigid-body modes can be written in the following two-partition form:  
+
+$$
+\boldsymbol{\Psi}_{r}\equiv\left[\frac{\boldsymbol{\Psi}_{i r}}{\boldsymbol{\Psi}_{e r}}\right]\equiv\left[\begin{array}{l}{\boldsymbol{\Psi}_{i r}}\\ {\boldsymbol{\Psi}_{b r}}\end{array}\right]
+$$  
+
+and then these rigid-body modes combined with the kept free-interface normal modes to form the following superset of kept modes:  
+
+$$
+\begin{array}{r}{\hat{\Phi}_{k}\quad\equiv\left[\begin{array}{l}{\hat{\Phi}_{i k}}\\ {\hat{\Phi}_{b k}}\end{array}\right]=\left[\begin{array}{l l}{\Psi_{i r}}&{\Phi_{i k}}\\ {\Psi_{b r}}&{\Phi_{b k}}\end{array}\right]}\end{array}
+$$  
+
+$$
+u^{c}\equiv\left\{\begin{array}{l}{{{\pmb u}_{i}}}\\ {{{\pmb u}_{b}}}\end{array}\right\}^{c}=\left[\begin{array}{l l}{{\hat{\Phi}}_{i k}}&{{\Psi_{i b}}}\\ {{\hat{\Phi}}_{b k}}&{{\Psi_{b b}}}\end{array}\right]^{c}\left\{\begin{array}{l}{{p_{k}}}\\ {{p_{b}}}\end{array}\right\}^{c}
+$$  
+
+where the “Rubin transformation matrir" is+  
+
+$$
+\begin{array}{r l r}{\Psi_{R}^{c}}&{=}&{\left[\begin{array}{l l}{\hat{\Psi}_{i r}}&{\hat{\Phi}_{i k}}\\ {\Psi_{b r}}&{\hat{\Phi}_{b k}}\end{array}\right|\,\Psi_{b b}\,\right]^{c}}\\ &{\equiv}&{\left[\begin{array}{l l}{\hat{\Phi}_{i k}}&{\Psi_{i b}}\\ {\hat{\Phi}_{b k}}&{\Psi_{b b}}\end{array}\right]^{c}}\end{array}
+$$  
+
+When two of the six free-interface modes for the 8DOF beam component are kept, the resulting eight columns of the Rubin transformation matrix have the shapes shown in Fig. 9 and 10.  
+
+![](36672dcc825d0016a59389a80204b6613410d02888b4831de1b6d08f762ca2f7.jpg)  
+Figure 9: Rubin Transformation Matrix-Cols. 1-4.  
+
+![](357f1a9c11f564c4cc1b03bf10872878921fb03e7dd2ab6b4a87004f5074a79d.jpg)  
+Figure 10: Rubin Transformation Matrix-Cols. 5-8.  
+
+With   mass-normalized   free-interface  normal modes $\Phi_{k}$ , from Eqs. (4) the component generalized mass matrix and stiffness matrix based on the Rubin transformation are  
+
+$$
+\begin{array}{r l}&{\bar{M}_{R}^{c}=\left[\begin{array}{l l l}{\bar{M}_{r r}}&{0_{r k}}&{0_{r b}}\\ {0_{k r}}&{I_{k k}}&{0_{k b}}\\ {0_{b r}}&{0_{b k}}&{\bar{M}_{b b}}\end{array}\right]^{c}}\\ &{\bar{K}_{R}^{c}\;=\left[\begin{array}{l l l}{0_{r r}}&{0_{r k}}&{0_{r b}}\\ {0_{k r}}&{\Lambda_{k k}}&{0_{k b}}\\ {0_{b r}}&{0_{b k}}&{\Psi_{b b}}\end{array}\right]^{c}}\end{array}
+$$  
+
+The zeros in these matrices are the result of orthogonality (e.g., Eqs. (37) and (38)), and the residual flexibility term $\Psi_{b b}$ of $\bar{K}_{R}^{c}$ is due to Eq. (35).s The above formulation is a consistent Ritz transformation; residual effects are included in both the stiffness and mass matrices. If the residual term $\bar{M}_{b b}$ in the mass matrix is omitted, the resulting nonconsistent transformation is referred to as the MacNeal method[6].  
+
+As suggested by Martinez, et al., let the lower partition of Eq. (65) be solved for the generalized coordinate vector $\pmb{p_{b}}$ in terms of $\mathbf{\deltau}_{b}$ , and the result incorporated back into the upper partition of Eq. (65). Finally, Eq. (65) can be re-cast in terms of the modal generalized coordinate vector $\pmb{p}_{k}$ and the interface displacement vector $\pmb{u}_{b}$ . This produces the following coordinate-transformation equation:  
+
+$$
+\begin{array}{r l r}{\lefteqn{\boldsymbol{u}^{c}\equiv\left\{\begin{array}{l}{\boldsymbol{u}_{i}}\\ {\boldsymbol{u}_{b}}\end{array}\right\}^{c}}}\\ &{=}&{\left[\begin{array}{l l}{\tilde{\Phi}_{i k}}&{\tilde{\Psi}_{i b}}\\ {0_{b k}}&{I_{b b}}\end{array}\right]^{c}\left\{\begin{array}{l}{\boldsymbol{p}_{k}}\\ {\boldsymbol{u}_{b}}\end{array}\right\}^{c}}\end{array}
+$$  
+
+This “Rubin-Martinez transformation matrir," has the same form as the Craig-Bampton transformation matrix of Eq. (58), and is given by  
+
+$$
+\begin{array}{r}{\everymath{\displaystyle}\Psi_{R M}^{c}=\left[\begin{array}{c c}{\tilde{\Phi}_{i k}}&{\tilde{\Psi}_{i b}}\\ {0_{b k}}&{I_{b b}}\end{array}\right]^{c}}\\ {=\left[\begin{array}{c c}{\left[\hat{\Phi}_{i k}-\Psi_{i b}\Psi_{b b}^{-1}\hat{\Phi}_{b k}\right]}&{\Psi_{i b}\Psi_{b b}^{-1}}\\ {0_{b k}}&{I_{b b}}\end{array}\right]^{c}}\end{array}
+$$  
+
+When two of the six free-interface modes for the 8DOF beam component are kept, the resulting eight columns of the Rubin-Martinez transformation matrix have the shapes shown in Figs. 11 and 12. Note the similarities and differences between the shapes in Figs. 11 and 3, and the similarities and differences between the shapes in Figs. 12 and 5.  
+
+When  the Rubin-Martinez  transformation  is used in Eqs. (4) to form the generalized mass matrix and stiffness matrix, the resulting matrices no longer have the sparsity that is exhibited by the Rubin mass matrix and stiffness matrix in Eqs. (67).  
+
+![](cb270e76b8f7c7d4b3d0875352e703251bbbedd3c3e1fb7cb0d0a7df5b6b5a63.jpg)  
+Figure 11: R-M Transformation Matrix-Cols. 1-4.  
+
+![](aab010052d671454c7da808acc33efd9e5de2c5009f3ebb9ab22b78f45832424.jpg)  
+Figure 12: R-M Transformation Matrix-Cols. 5--8.  
+
+# 5. Conclusions and Recommendations  
+
+Procedures used to formulate component modes for substructures and to assemble substructure models to form reduced-order models of the original system have been reviewed. The physical meaning of many CMS terms has been illustrated. The constraint-mode method described in Section 4.1 is widely used in reducing finite element models for dynamic analysis because it is very straightforward, and it leads to accurate reduced-order models.  On the other hand, the attachment-mode method described in Section 4.2 also produces accurate reduced-order models, and is currently widely used in test-verifying finite element models. There is a pressing need for the development of efficient computational structural dynamics algorithms based on constraint-mode substructuring methods, and there is still a substantial need for a more thorough understanding of attachment-mode methods and their use in test verification of finite element models. Research is also needed on CMS methods for damped structural systems.  
+
+# 6. References  
+
+[1] Bisplinghoff, R. L., Ashley, H., and Halfman, R. L., Aeroelasticity, Addison-Wesley Publishing Co., Inc., Reading, Ma, 1955.   
+[2] Hurty, W. C., Dynamic Analysis of Structural Systems by Component Mode Synthesis, Tech. Report No. 32-530, Jet Propulsion Laboratory, Pasadena, CA, Jan. 1964.   
+[3] Hurty, W. C., “Dynamic Analysis of Structural Systems Using Component Modes," AIAA Journal, Vol. 3, No. 4, 1965, pp. 678-685.   
+[4] Bamford, R. M., A Modal Combination Program for Dynamic Analysis of Structures (Revision No. 1), Tech. Memo. No. 33-290, Jet Propulsion Laboratory, Pasadena, CA, July 1, 1967.   
+[5] Craig, R. R., Jr., and Bampton, M. C. C., "Coupling of Substructures for Dynamic Analysis," AIAA Journal, Vol. 6, No. 7, 1968, Pp. 1313-1319.   
+[6] MacNeal, R. H., “"A Hybrid Method of Component Mode Synthesis," J. Computers & Structures, Vol. 1, No. 4, Dec. 1971, pp. 581-601.   
+[7] Rubin, S., “Improved Component-Mode Representation for Structural Dynamic Analysis," AIAA Journal, Vol. 13, No. 8, Aug. 1975, pp. 995-1006.   
+[8] Benfield, W. A., Bodley, C. S., and Morosow, G., “Modal Synthesis Methods," Space Shuttle Dynamics and Aeroelasticity Working Group Symposium on Substructure Testing and Synthesis, NASA-TM-X-72318, Marshall Space Flight Center, AL, 1972.  
+
+[9] Craig, R. R., Jr., "A Review of Time-Domain and Frequency-Domain Component Mode Synthesis Methods,” Int. J. Analytical and Ecperimental Modal Analysis, Vol. 2, No. 2, 1987, pp.59-72.  
+
+[10] Craig, R. R., Jr., “Substructure Methods in Vibration," ASME Transactions, Special 50th Anniversary Design Issue, Vol. li7, June 1995, pp. 207-213.   
+[11] Craig, R. R., Jr., Structural Dynamics - An Introduction to Computer Methods, John Wiley & Sons, Inc., New York, NY, 1981.   
+[12] Meirovitch, L., Computational Methods in Structural Dynamics, Sijthoff & Noordhoff, Rockville, MD, 1980.   
+[13] Maia, N. M. M., and Silva, J. M. M., eds., Theoretical and Eaperimental Modal Analysis,(Research Studies Press, Ltd.), John Wiley & Sons, Inc., New York, 1997.   
+[14] Martinez, D. R., et al., “Combined Experimental/Analytical Modeling Using Component Mode Synthesis," AIAA Paper 84- 0941, AIAA/ASME/ASCE/AHS 25th Structures, Structural Dynamics and Materials Conf, Palm Springs, CA, May 1984, pp. 140- 152.   
+[15] Kammer, D. C., and Baker, M., “A Comparison of the Craig-Bampton and Residual Flexibility Methods for Component Substructure Representation," AIAA Journal of Aircraft, Vol. 24, No. 4, March 1987, pp. 262-267.   
+[16] Admire, J. R., Tinker, M. L., and Ivey, E. W., "Mass-Additive Modal Test Method for Verification of Constrained Structural Models," AIAA Journal, Vol. 31, No. 11, 1993, pp. 2148- 2153.   
+[17]  Admire, J. R., Tinker, M. L., and Ivey, E. W., "Residual Flexibility Test Methods for Verification of Constrained Structural Models," AIAA Journal, Vol. 32, No. 1, Jan. 1994, pp. 170-175.   
+[18] Duarte, M. L. M., and Ewins, D. J., “Improved Experimental Component Mode Synthesis (IECMS) with Residual Compensation Based Purely on Experimental Results," Proc. 14th International Modal Analysis Conference, Dearborn, MI, Feb. 1996, pp. 641-647.   
+[19] Chandler, K. O., and Tinker, M. L., “A General MassAdditive Method for Component Mode Synthesis," Paper No. AIAA-97-1381, Proc. 38th Structures, Structural Dynamics and Materials Conf., Kissimmee, FL, Apr. 1997, pp. 93- 103.   
+[20] Tinker, M. L., and Cutchins, M. A., “Model Correlation Issues in Residual Flexibility Testing," Paper No. DETC97/VIB-4262, Proc. 1997 ASME Design Engineering Technical Conferences, Sacramento, CA, Sept. 1997.   
+[21] Bennighof, J. K., Kaplan, M. F., and Muller, M. B., “Extending the Frequency Response Capabilities of Automated Multi-Level Substructuring," AIAA Dynamics Specialists Conference, Atlanta, GA, April 2000, to appear.   
+[22] Farhat, C., Lesoinne, M., and Pierson, K., “A Scalable Substructuring Method for Transient and Vibration Analyses on Massively Parallel Processors," AIAA Dynamics Specialists Conference, Atlanta, GA, April 2000, to appear.   
+[23] Craig, R. R., Jr., and Hale, A. L., “BlockKrylov Component Synthesis Method for Structural Model Reduction,” AIAA J. Guidance, Control, and Dynamics, Vol. 11, No. 6, 1988, pp. 562-570.   
+[24] Benfield, W. A., and Hruda, R. F., “Vibration Analysis of Structures by Component Mode Substitution,” AIAA Journal, vol. 9, No. 7, July 1971, pp. 1255-1261.   
+[25] Hintz, R. M., “Analytical Methods in Component Modal Synthesis,” AIAA Journal, Vol. 13, No. 8, 1975, pp. 1007-1016.   
+[26] Craig, R. R., Jr., and Chang, C-J., "On the Use of Attachment Modes in Substructure Coupling for Dynamic Analysis," Paper 77-405, Proc. 18th Structures, Structural Dynamics and Materials Conf., San Diego, CA, 1977, pp. 89-99.   
+[27] “Introduction to Superelements in Dynamic Analysis," MSC/NASTRAN REFERENCE MANUAL, Ver. 69, Chapt. 10, The MacNealSchwendler Corp. Los Angeles, CA, 1996.   
+[28] Martinez, D. R., and Gregory, D. L., A Comparison of Free Component Mode Synthesis Techniques Using MSC/NASTRAN, Report No. SAND83-0025, Sandia National Laboratories, Albuquerque, NM, June 1984.  

多体+耦合求解器/参考文献/Nim-Section force calculation in flexible subs/auto/Nim-Section force calculation in flexible subs.md

@@ -0,0 +1,265 @@
+PAPER • OPEN ACCESS  
+
+# Section force calculation in flexible substructures modelled by wind turbine design tool Bladed  
+
+To cite this article: Erik Nim et al 2024 J. Phys.: Conf. Ser. 2767 052043  
+
+View the article online for updates and enhancements.  
+
+You may also like  
+
+-Correction of periodic displacement nonlinearities by two-wavelength   
+interferometry   
+Angus Bridges, Andrew Yacoot, Thomas Kissinger et al.   
+Odd Random Phase Electrochemical   
+Impedance Spectroscopy to Study the Corrosion Behavior of Hot Dip Zn and ZnAlloy Coated Steel Wires in Sodium   
+Chloride Solution   
+Gopal Ji, Lucía Fernández Macía, Bart Allaert et al.   
+-New evidence and impact of electron   
+transport non-linearities based on new perturbative inter-modulation analysis M. van Berkel, T. Kobayashi, H. Igami et al.  
+
+![](fc3dbed62a6f9f12efeebd7537e4aab4a605f39df121cf893f293d0b3b6a8f8a.jpg)  
+
+# Section force calculation in flexible substructures modelled by wind turbine design tool Bladed  
+
+Erik $\mathbf{Nim}^{1}$ , John Roadnight2, Hassan Moharram2, William Collier2 and Chr. Sigurd L. Jensen  
+
+1 DNV Denmark A/S, Ecopark, Bautavej 1A, 8210 Aarhus V, Denmark   
+2 DNV Services UK Limited, One Linear Park, Bristol, BS2 0PS, United Kingdom  
+
+E-mail: erik.nim@dnv.com  
+
+Abstract. Aeroelastic simulation tools are widely used to predict the coupled response for onshore, offshore and floating wind turbines. A key output of the analysis is the section forces in the flexible substructures, such as blades, towers, and floating foundations. Individual substructures are typically modelled as a single linear flexible body, which is a part of a multibody system used for modelling the complete wind turbine. However, in some floating foundation concepts, in particular those with pre-tensioned cables, geometric non-linearities in the substructure can have a significant effect. With the conventional method for calculating section forces, such non-linearities are not captured accurately. In this paper, a new equilibriumbased method for calculation of section forces is presented, as part of the wind turbine design tool Bladed. This method relies on equilibrium between the section forces and the applied loads. The section forces are determined in terms of Lagrange multipliers, which represent the internal constraint forces in the structural elements. This method is valid even when significant geometric non-linearities are present, with reduced order models, and for statically indeterminate substructures. A case study is presented demonstrating the self-consistency of the new method for the OC4 semi-submersible floating platform, modified to include pre-tensioned cables. A verification against ANSYS is also presented for some elementary cases.  
+
+# 1. Introduction and motivation  
+
+Aeroelastic simulation tools are widely used to predict the coupled response for onshore, offshore and floating wind turbines. In wind turbine design tool Bladed [1], the structural model is based on a multibody system approach, where the complete wind turbine structure is modelled as an assembly of rigid and flexible bodies representing typical turbine substructures such as blades and tubular towers, as well as complex offshore support structures such as jackets and floating foundations.  
+
+A key functionality for a flexible body is the calculation of interior section forces (also known as internal forces). In the conventional method for calculating section forces, the calculation is carried out in the undeflected configuration. This is a good assumption in many cases but can lead to inaccuracies in bending and torsion moments when geometrical non-linearities become more significant. The inaccuracy is observed most clearly when calculating section forces in floating foundations, which can have significant non-linearities introduced by certain features, e.g., pre-tensioned cables. The present work represents a novel contribution to address these inaccuracies in an efficient and consistent way.  
+
+This paper proposes a new equilibrium-based method for the calculation of section forces in flexible bodies. The theoretical approach is described in section 2. The implementation in Bladed is verified against ANSYS in section 3, and a case study is presented for a floating substructure with pre-tensioned cables in section 4.  
+
+# 2. Theory summary  
+
+In this section, a new method for calculating section forces in flexible bodies is described. The method is based on an equilibrium approach often used in ad-hoc manual analysis of simple problems, but it has also been used as a basis for more general methods applicable to both statically determinate and indeterminate space frames [2][3].  
+
+# 2.1. Modelling of flexible bodies  
+
+In Bladed, the flexible bodies are modelled as structural frames using a linear finite element method [4][5] combined with the multibody system theory described by Shabana [6]. A structural frame is modelled by conventional finite element technique using two-node Timoshenko beam elements with 12 degrees of freedom (DOFs) along with a few additional elements for modelling specific structural parts, such as rigid bodies and bar elements used for modelling cables. An explicit geometric stiffness model is introduced to account for the second-order effects of the applied loads [4], but this approach has some limitations and is generally not sufficiently accurate. A traditional Craig-Bampton type modal truncation method is introduced to reduce the total number of degrees of freedom of the flexible body and to filter out unwanted high-frequency modes [7][8].  
+
+It is well known that the introduction of a modal approach for system reduction in multibody systems introduces a problem for representing the section forces sufficiently accurately, and it has therefore been proposed to calculate the section forces in a post-processing step, e.g., using finite element codes [9][10]. In Bladed, the section forces in the flexible bodies are calculated in a separate step at the modal truncated state resulting from dynamic or static analyses. For dynamic analyses, this post-processing step must be performed after every time step, which means that computational performance is critical.  
+
+# 2.2. Conventional calculation of section forces with the finite element method  
+
+With the displacement-based finite element method, the section forces are calculated at the end nodes of the beam elements using the element stiffness matrix. Possible geometric constraints, specified by the user to constrain some displacement patterns such as the axial extension of tubular towers, are modelled by a set of linear constraint equations with an associated set of Lagrange multipliers.  
+
+The conventional linear method for calculating section forces does not include any second-order effects of the loads applied to the flexible body due to fact that the loads are applied to the reference state and does not account for any deflections. This introduces an inaccuracy in the calculated bending and torsion moments, which can introduce significant errors in some load components, e.g., the torsion moment along a highly flexible blade. A geometric stiffness matrix is added to the material stiffness matrix to include some second-order effects [4].  
+
+# 2.3. A new equilibrium-based method for section force calculation  
+
+The most simple and accurate method for the calculation of section forces in a flexible body at a known deformed state is derived by equilibrating the section forces with the applied loads. This method allows the calculation to be carried out using only the loading, the actual deflection state and the geometry as well as the element connectivity, without using any stiffness (or damping) properties. However, the method only applies to statically determinate flexible bodies such as blades and tubular towers. For indeterminate flexible bodies including most support structures, the section forces depend on the stiffness properties, which consequently must be included in the calculation.  
+
+The new equilibrium-based method employs the multibody system approach used by Bladed for modelling of the complete wind turbine structure [1][11]. With this approach, a flexible body is modelled as an assembly of rigid and flexible elements that are interconnected at a number of userdefined nodes. A subset of these nodes is selected as global nodes used for connecting the flexible body to the neighbouring structure, whereas the remaining interior nodes are only used internally.  
+
+In the following, it is assumed that the deformation of a flexible body element is described by a set of generalised strains, which for a two-node beam element without constraints is conveniently collected in a vector $\pmb{\varepsilon}^{\mathrm{e}}$ . In case that any deformation mode, such as torsion or axial extension, is constrained, the number of generalised strains is reduced accordingly. The absolute position and orientation of a node is described by a position vector and a rotation matrix to enable description of finite rotations. It is important to note that these quantities are dependent variables that are determined by the generalised strains. Although it is not possible to describe the absolute position and orientation of the element nodes by a vector because finite rotations are not vector quantities, a variation, i.e., a small virtual change, of these quantities can be described by a common vector $\delta\mathbf{x}_{\mathrm{r}}^{\mathrm{e}}$ , which for a two-node beam element is a 12-by-1 vector. The relation between a variation $\delta\mathbf{x}_{\mathrm{r}}^{\mathrm{e}}$ of the nodal position and orientation and a variation $\delta\pmb{\varepsilon}^{\mathrm{e}}$ of the generalised strains defines the geometric constraints for the element. These constraint relations are written in a linear form  
+
+$$
+\mathbf{C}_{\mathrm{r}}^{\mathrm{e}}\,\delta\mathbf{x}_{\mathrm{r}}^{\mathrm{e}}+\mathbf{C}_{\mathrm{\varepsilon}}^{\mathrm{e}}\,\delta\mathbf{\varepsilon}^{\mathrm{e}}=\mathbf{o}
+$$  
+
+where $\mathbf{C}_{\mathrm{r}}^{\mathrm{e}}$ and $\mathbf{C}_{\varepsilon}^{\mathrm{e}}$ are geometric constraint matrices that generally depend on the deformation of the element. The symbol 𝐨 introduced in (1) represents a null vector with matching dimension. In general, the constraint matrices must be derived for all types of elements based on the geometry and the internal kinematics of the element. The derivation of the expressions for the constraint matrices for the used beam element was done using the co-rotational approach based on six natural deformation modes [12].  
+
+The derivation of the dynamic equilibrium equations for a flexible body element is based on the principle of virtual work similar to what is done for the multibody system modelling the complete structure [1][11]. It turns out that the element contribution to the virtual work can be written in the simplified form  
+
+$$
+\delta W^{\mathrm{e}}=\left[\begin{array}{l}{\delta\mathbf{x}_{\mathrm{r}}^{\mathrm{e}}}\\ {\delta\mathbf{\varepsilon}^{\mathrm{e}}}\end{array}\right]^{T}\left(\left[\begin{array}{l}{\mathbf{C}_{\mathrm{r}}^{\mathrm{e}^{T}}}\\ {\mathbf{C}_{\varepsilon}^{\mathrm{e}T}}\end{array}\right]\lambda^{\mathrm{e}}-\left[\begin{array}{c}{\mathbf{p}_{\mathrm{r}}^{\mathrm{e}}+\mathbf{f}_{0}^{\mathrm{e}}}\\ {\mathbf{p}_{\varepsilon}^{\mathrm{e}}-\mathbf{K}_{\varepsilon\varepsilon}^{\mathrm{e}}\,\varepsilon^{\mathrm{e}}}\end{array}\right]\right)
+$$  
+
+where the vectors ${\bf p}_{\mathrm{r}}^{\mathrm{e}}$ and ${\bf p}_{\varepsilon}^{\mathrm{e}}$ represent the applied element loads including inertial loads conjugate to respectively nodal motion and generalised strain, while the stiffness matrix ${\bf K}_{\varepsilon\varepsilon}^{\mathrm{e}}$ represent the elastic properties on the assumption of linear elastic material. In addition, the vector $\lambda^{\mathrm{e}}$ represents the yet unknown Lagrange multipliers associated with the constraint matrices, while the vector $\mathbf{f}_{0}^{\mathrm{e}}$ contains the element section forces, which originate from the connection to neighbouring elements. The element equilibrium equations corresponding to the virtual work expression (2) is found by requiring $\delta W^{\mathrm{e}}=0$ for any combination of the variations $\delta\mathbf{x}_{\mathrm{r}}^{\mathrm{e}}$ and $\delta\pmb{\varepsilon}^{\mathrm{e}}$ :  
+
+$$
+\left[\begin{array}{c}{\mathbf{C}_{\mathrm{r}}^{\mathrm{e}^{T}}}\\ {\mathbf{C}_{\mathrm{\varepsilon}}^{\mathrm{e}T}}\end{array}\right]\boldsymbol{\lambda}^{\mathrm{e}}=\left[\begin{array}{c}{\mathbf{p}_{\mathrm{r}}^{\mathrm{e}}+\mathbf{f}_{0}^{\mathrm{e}}}\\ {\mathbf{p}_{\mathrm{\varepsilon}}^{\mathrm{e}}-\mathbf{K}_{\mathrm{\varepsilon}\mathrm{\varepsilon}}^{\mathrm{e}}\varepsilon^{\mathrm{e}}}\end{array}\right]
+$$  
+
+The resulting system for the complete flexible body is derived with the conventional finite element technique by assembling the element equations. To this end, the variation of the nodal positions and orientations for the complete flexible body are collected in a common $6N.$ -by-1 vector $\delta\mathbf{x}_{\mathrm{r}}$ , where $N$ denotes the total number of nodes of the flexible body. In addition, all generalised strain components are collected in a common $N_{\varepsilon}$ -by-1 strain vector $\pmb{\varepsilon}$ , where $N_{\varepsilon}$ denotes the total number of generalised strains, while the Lagrange multipliers are collected in a common $N_{\mathsf{c}}$ -by-1 vector $\lambda$ , where $N_{\mathsf{c}}$ denotes the total number of geometric constraints.  
+
+The geometric constraints for the complete flexible body are found by assembling the element constraints (1):  
+
+$$
+\mathbf{C}_{\mathrm{r}}\,\delta\mathbf{x}_{\mathrm{r}}+\mathbf{C}_{\varepsilon}\,\delta\mathbf{\varepsilon}=\mathbf{o}
+$$  
+
+where $\mathbf{C}_{\mathrm{r}}$ is the $N_{\mathsf{c}}$ -by- $.6N$ constraint matrix for nodal displacement and rotation, while $\mathbf{C}_{\varepsilon}$ is the $N_{\mathsf{c}}$ -by$N_{\varepsilon}$ constraint matrix for the generalised strains.  
+
+The virtual work for the complete flexible body simply becomes $\delta W=\textstyle\sum_{\mathrm{e}}\delta W^{\mathrm{e}}$ , where $\delta W^{\mathbf{e}}$ are the element virtual work contribution (2). It is straightforward to show that the resulting virtual work contribution can be written in the form  
+
+$$
+\delta W=\left[\begin{array}{c}{\delta\mathbf{x}_{\mathrm{r}}}\\ {\delta\mathbf{\varepsilon}}\end{array}\right]^{T}\left(\left[\begin{array}{c}{\mathbf{C}_{\mathrm{r}}^{\mathbf{\Gamma}}}\\ {\mathbf{C}_{\varepsilon}^{\mathbf{\Gamma}}\mathbf{\Gamma}}\end{array}\right]\lambda\mathbf{\varepsilon}-\left[\begin{array}{c}{\mathbf{p}_{\mathrm{r}}}\\ {\mathbf{p}_{\varepsilon}-\mathbf{K}_{\varepsilon\varepsilon}\mathbf{\varepsilon}}\end{array}\right]\right)
+$$  
+
+where the flexible body load vectors $\mathbf{p_{r}}$ and $\mathbf{p_{\varepsilon}}$ , stiffness matrix $\mathbf{K}_{\varepsilon\varepsilon}$ , and Lagrange multipliers $\lambda$ are found by assembling the corresponding element quantities. The flexible body equations corresponding to the virtual work contribution (5) becomes  
+
+$$
+\left[\begin{array}{c}{\mathbf{C}_{\mathrm{r}}^{\mathrm{\Delta}T}}\\ {\mathbf{C}_{\mathrm{\varepsilon}}^{\mathrm{\Delta}T}}\end{array}\right]\boldsymbol{\lambda}=\left[\begin{array}{c}{\mathbf{p}_{\mathrm{r}}}\\ {\mathbf{p}_{\mathrm{\varepsilon}}-\mathbf{K}_{\mathrm{\varepsilon}\varepsilon}\;\varepsilon}\end{array}\right]
+$$  
+
+# 2.3.1. Statically determinate flexible bodies  
+
+For statically determinate flexible bodies, the number of constraints equals the number of nodal DOFs, i.e. 6 times the number of nodes. In this case, the Lagrange multipliers $\lambda$ for the flexible body are calculated from the linear system $\mathbf{C}_{\mathrm{r}}^{\textit{T}}\hat{\mathbf{\alpha}}=\mathbf{p}_{\mathrm{r}}$ appearing from (6) with the assumption that the constraint matrix $\mathbf{C}_{\mathrm{r}}$ is invertible. The extracted element Lagrange multipliers $\lambda^{\mathrm{e}}$ are then used for calculating the resulting section forces at the element end nodes by the linear system ${\bf f}_{0}^{\mathrm{e}}={\bf C}_{\mathrm{r}}^{\mathrm{e}^{T}}{\bf\lambda}^{\mathrm{e}}-{\bf p}_{\mathrm{r}}^{\mathrm{e}}$ derived from (3). With this approach, the section forces are determined in terms of the Lagrange multipliers, which represent the internal constraint forces in the elements.  
+
+It is important to note that the flexible body constraint matrix $\mathbf{C_{\mathrm{r}}}$ in (4) and the element constraint matrix $\mathbf{C}_{\mathrm{r}}^{\mathrm{e}}$ in (1) must be evaluated for the modal truncated deflection state used during time integration or the steady state analysis of the complete system. In particular, the modal truncation can exclude all modes in which case the flexible body appears to be rigid.  
+
+# 2.3.2. Statically indeterminate flexible bodies  
+
+For a statically indeterminate flexible body, the number of constraints is larger than the number of nodal coordinates, i.e., $N_{\mathrm{c}}>6N$ , which means that it is not possible to calculate the section forces from equilibrium alone, as is the case for a statically determinate flexible body. The degree of static indeterminacy, in the following denoted as $D_{s}$ , is commonly used for quantifying the number of redundant forces, which in the present context means that $D_{\mathrm{s}}=N_{\mathrm{c}}-6N$ . It is clear that ${{D}_{s}}=0$ for a statically determinate flexible body, while $D_{s}>0$ in the statically indeterminate case.  
+
+The approach for calculating section forces in statically indeterminate flexible bodies is based on subdividing the $N_{\mathsf{c}}$ geometric constraints for a flexible body into a set of $6N$ determinate constraints and the remaining $D_{\mathrm{s}}$ indeterminate constraints. In the following, the Lagrange multipliers associated with the determinate constraints are collected in the 6N-by-1 vector $\lambda_{\mathrm{r}}$ , while the Lagrange multipliers associated with the remaining indeterminate constraints are collected in the $D_{\mathrm{s}}$ -by-1 vector $\lambda_{\mathrm{e}}$ . The subdivision of the Lagrange multipliers and hence the constraint relations is conveniently described in terms of a permutation matrix $[\mathbf{P}_{\mathrm{r}}\;\mathbf{P}_{\mathrm{e}}]$ in the form  
+
+$$
+\lambda=\left[\mathrm{\bf~P}_{\mathrm{r}}\;\mathrm{\bf~P}_{\mathrm{e}}\right]\left[\lambda_{\mathrm{r}}\right]=\mathrm{\bf~P}_{\mathrm{r}}\,\lambda_{\mathrm{r}}+\mathrm{\bf~P}_{\mathrm{e}}\,\lambda_{\mathrm{e}}
+$$  
+
+where $\mathbf{P}_{\mathrm{r}}$ is a $N_{\mathsf{c}}$ -by- $.6N$ partial permutation matrix for identifying the determinate constraints while $\mathbf{P_{e}}$ is a $N_{\mathsf{c}}$ -by- $\cdot D_{s}$ partial permutation matrix for identifying the indeterminate constraints. The corresponding subdivision of the constraint relations (4) results in the system  
+
+$$
+\begin{array}{r}{\mathbf{C}_{\mathrm{rr}}\,\delta\mathbf{x}_{\mathrm{r}}+\mathbf{C}_{\mathrm{r}\varepsilon}\,\delta\mathbf{\varepsilon}=\mathbf{o}}\\ {\mathbf{C}_{\mathrm{er}}\,\delta\mathbf{x}_{\mathrm{r}}+\mathbf{C}_{\mathrm{e}\varepsilon}\,\delta\mathbf{\varepsilon}=\mathbf{o}}\end{array}
+$$  
+
+where $\mathbf{C}_{\mathrm{rr}}=\mathbf{P}_{\mathrm{r}}^{\,\,T}\mathbf{C}_{\mathrm{r}}$ and ${\bf C}_{\mathrm{r}\mathrm{s}}={\bf P}_{\mathrm{r}}^{~T}{\bf C}_{\mathrm{s}}$ , while $\mathbf{C}_{\mathrm{er}}=\mathbf{P}_{\mathrm{e}}^{\ \ T}\mathbf{C}_{\mathrm{r}}$ and ${\bf C}_{\bf e\mathrm{s}}={\bf P}_{\bf e}^{\mathrm{~\tiny~T~}}{\bf C}_{\bf s}$ . With this subdivision, the virtual work expression (5) takes the form  
+
+$$
+\delta W=\left[\begin{array}{c}{\delta\mathbf{x_{r}}}\\ {\delta\varepsilon}\end{array}\right]^{T}\left(\left[\begin{array}{c c}{\mathbf{C_{rr}}^{T}}&{\mathbf{C_{er}}^{T}}\\ {\mathbf{C_{r\varepsilon}}^{T}}&{\mathbf{C_{e\varepsilon}}^{T}}\end{array}\right]\left[\begin{array}{c}{\lambda_{\mathrm{r}}}\\ {\lambda_{\mathrm{e}}}\end{array}\right]-\left[\begin{array}{c}{\mathbf{p_{r}}}\\ {\mathbf{p_{\varepsilon}}-\mathbf{K_{\varepsilon\varepsilon}}\,\varepsilon}\end{array}\right]\right)
+$$  
+
+In principle, the subdivision of the geometric constraints into determinate and indeterminate partitions corresponds to the introduction of a set of virtual structural cuts that makes an indeterminate flexible body statically determinate and enables the calculation of section forces by an equilibrium approach. It is assumed that the subdivision of the constraints is carried out in such a way that the square matrix $\mathbf{C}_{\mathrm{rr}}$ is invertible, which ensures that the constraint equation (8) can be solved with respect to the variation $\delta\mathbf{x}_{\mathrm{r}}$ of the nodal position and orientation. With this assumption, it follows that the subdivided constraint relations (8) and (9) can be written as  
+
+$$
+\delta\mathbf{x}_{\mathrm{r}}=\mathbf{T}_{\mathrm{r}\mathrm{{s}}}\delta\mathbf{\varepsilon}
+$$  
+
+$$
+\bigl(\mathbf{C}_{\mathrm{e}\varepsilon}+\mathbf{C}_{\mathrm{er}}\,\mathbf{T}_{\mathrm{r}\varepsilon}\bigr)\delta\mathbf{\varepsilon}=\mathbf{o}
+$$  
+
+where $\mathbf{T}_{\mathbf{r}\mathbf{s}}$ is a transformation matrix that is calculated from the linear system $\mathbf{C}_{\mathrm{rr}}\,\mathbf{T}_{\mathrm{r\varepsilon}}=-\mathbf{C}_{\mathrm{r}\mathrm{\varepsilon}}$ .  
+
+In order to eliminate the variation $\delta\mathbf{x}_{\mathrm{r}}$ of the nodal position and orientation in the virtual work expression (10), the transformation (11), together with the trivial identity $\delta\pmb{\varepsilon}=\delta\pmb{\varepsilon}$ , are written in the common linear form  
+
+$$
+\left[\begin{array}{l}{\delta\mathbf{x}_{\mathrm{r}}}\\ {\delta\mathbf{\varepsilon}}\end{array}\right]=\left[\begin{array}{l}{\mathbf{T}_{\mathrm{r}\varepsilon}}\\ {\mathbf{I}}\end{array}\right]\delta\mathbf{\varepsilon}
+$$  
+
+where the symbol 𝐈 represents the identity matrix with matching dimensions. With this transformation, the virtual work expression (10) can be written in terms of the generalised strain variations $\delta\pmb{\varepsilon}$ as  
+
+$$
+\delta W=\delta\varepsilon^{T}\left(\left(\mathbf{C_{e\varepsilon}}+\mathbf{C_{er}}\,\mathbf{T_{r\varepsilon}}\right)^{T}\lambda_{\mathrm{e}}-\left(\mathbf{p_{\varepsilon}}+\mathbf{T_{r\varepsilon}}^{T}\mathbf{p_{r}}-\mathbf{K_{\varepsilon\varepsilon}}\varepsilon\right)\right)
+$$  
+
+The resulting system of equilibrium equations corresponding to the virtual work expression (14) is found by requiring $\delta W=0$ for any combination of the generalised strain variations $\delta\pmb{\varepsilon}$ . It is convenient to write the resulting system together with the constraint relation (12) in terms of a linear system  
+
+$$
+\left[\begin{array}{c c}{\mathbf{K}_{\mathrm{{e}}\varepsilon}^{*}}&{\mathbf{C}_{\mathrm{e}\varepsilon}^{*\ {T}}}\\ {\mathbf{C}_{\mathrm{e}\varepsilon}^{*}}&{\mathbf{0}}\end{array}\right]\left[\begin{array}{l}{\pmb{\varepsilon}^{\prime}}\\ {\lambda_{\mathrm{e}}}\end{array}\right]=\left[\begin{array}{l}{\mathbf{p}_{\mathrm{{e}}}^{*}}\\ {\mathbf{o}}\end{array}\right]
+$$  
+
+where ${\bf C}_{\mathrm{e}\varepsilon}^{\ast}={\bf C}_{\mathrm{e}\varepsilon}+{\bf C}_{\mathrm{er}}{\bf T}_{\mathrm{r}\varepsilon}$ and $\mathbf{p}_{\mathrm{{z}}}^{\ast}=\mathbf{p}_{\mathrm{{z}}}+\mathbf{T}_{\mathrm{{r}{\boldsymbol{\mathrm{\varepsilon}}}}}^{\;\;\;T}\mathbf{p}_{\mathrm{{r}}}$ , while $\mathbf{K}_{\varepsilon\varepsilon}^{*}=\mathbf{K}_{\varepsilon\varepsilon}$ for notational consistency. This linear system enables the calculation of the indeterminate Lagrange multipliers $\lambda_{\mathrm{e}}$ as well as the generalised strains $\pmb\varepsilon^{\prime}$ corresponding to the applied loads ${\bf p}_{\varepsilon}^{*}$ . Here the symbol $\pmb\varepsilon^{\prime}$ is introduced to emphasise that the generalised strains calculated from the linear system (15) define a refined displacement field rather than the modal truncated displacement field used in the static or the dynamic analysis of the complete wind turbine. In particular, this approach allows an accurate calculation of the section forces for simulations with a low number of mode shapes (possibly none) for describing elastic deformations of the flexible body component.  
+
+In principle, the linear system (15) can be solved directly using Gaussian elimination, but it turns out that the augmented stiffness matrix (i.e., the coefficient matrix) on the left-hand side is ill-conditioned. A numerically stable solution method was derived based on the constraint elimination method [5], where the generalised strains $\pmb\varepsilon^{\prime}$ are subdivided into a set of strains to be eliminated and a set of remaining independent strains. With this approach, the linear system (15) is transformed into a well-conditioned system, where the only unknowns are the independent generalised strains.  
+
+The determinate Lagrange multipliers $\lambda_{\mathrm{r}}$ corresponding to the vector $\lambda_{\mathrm{e}}$ of indeterminate Lagrange multipliers is calculated from the linear system ${{\bf{C}}_{\mathrm{{rr}}}}^{T}{\lambda_{\mathrm{{r}}}}={{\bf{p}}_{\mathrm{{r}}}}-{{\bf{C}}_{\mathrm{{er}}}}^{T}{\lambda_{\mathrm{{e}}}}$ that is found from the virtual work expression (10) by requiring $\delta W=0$ for any combination of the variations $\delta\mathbf{x}_{\mathrm{r}}$ . The complete vector $\lambda$ of Lagrange multipliers can then be assembled by (7), which means that the section forces can be calculated similarly to what is done in the statically determinate case.  
+
+# 3. Verification against ANSYS  
+
+Figure 1 illustrates the general form of the elementary test models that are used for comparison between the section forces computed by Bladed and ANSYS. The models are statically indeterminate structural frames comprising of 4 cylindrical beam members interconnected at 4 vertices.  
+
+![](a74f98e89ecb3807c00a405a0e09aa01094ae9053a9aeceee8b63cb90f0da896.jpg)  
+Figure 1. General definition of elementary verification model. Note that S1 to S5 locations for section force extraction are coincident with member ends.  
+
+As indicated in Figure 1, two different model configurations (denoted as 1 and 2) have been considered. In Model 2, the right and left vertices (V2 and V4) are joined by the addition of a bar element with no rotational degrees of freedom at each end.  
+
+All members are assigned tubular sectional properties with the outer diameter of $4\;\mathrm{m}$ and $20\;\mathrm{mm}$ thickness, with steel material with elastic modulus $E=200$ GPa, Poisson’s ratio $\mathbf{0}=0.3$ , shear modulus $G=76.9\;\mathrm{GPa}$ .  
+
+Three different load cases are solved in ANSYS and Bladed, as listed in Table 1, representing different combinations of concentrated point loads and distributed line loads applied to the models.  
+
+The Bladed load cases are individually solved as static non-linear analyses based on the Bladed multibody system in a separate application. All load cases in ANSYS are individually solved as static non-linear analyses, with loads applied over a single load step.  
+
+Bladed is solved with the geometric non-linearity option of ‘Support Structure Geometric Stiffness Model’ set as ‘Axial loads only’ [1]. Geometric non-linearity is switched on in ANSYS using the ‘NLGEOM’ option [13].  
+
+Nodal deflections are extracted at each of the vertices V1 to V4. Section forces and moments are extracted from the member section locations labelled S1 to S5 in Figure 1. These section locations are located immediately adjacent to the vertices.  
+
+Table 1. Definition of load cases. All loads are given in the global frame   
+
+
+<html><body><table><tr><td rowspan="2">Model</td><td rowspan="2">Load case</td><td colspan="3">Point loads on V2 (kN) (Right vertex)</td><td colspan="3">Point loads on V3 (kN) (Top vertex)</td><td colspan="3">Line loads on Mbr2 (kN/m) (Upper right member)</td><td colspan="3">Line loads on Mbr3 (kN/m) (Upper left member)</td></tr><tr><td>Px</td><td>Py</td><td>Pz</td><td>Px</td><td>Py</td><td>Pz</td><td>Lx</td><td>Ly</td><td>Lz</td><td>Lx</td><td>Ly</td><td>Lz</td></tr><tr><td>1</td><td>1A</td><td>50,000</td><td>5000</td><td>500</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td></td><td>0</td></tr><tr><td>1</td><td>1B</td><td>0</td><td>0</td><td>0</td><td>100,000</td><td>10,000</td><td>0</td><td>1000</td><td>1</td><td>1</td><td>-1000</td><td>0 -1</td><td>-1</td></tr><tr><td>2</td><td>2A</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>-1000</td><td>0</td><td>-1</td><td>-1000</td><td>0</td><td>1</td></tr></table></body></html>  
+
+# 3.1. Beam element discretisation  
+
+Each frame member of the verification models is defined in Bladed using a single two-node beam element. The outer frame members in the ANSYS models are meshed with three-node BEAM189 elements [13] with two element divisions per frame member. For Model 2, there is a horizontal bar element between the left and right vertices. This bar is modelled using a LINK180 element in ANSYS [13]. This element has three translational degrees of freedom at each end node, and no rotational stiffness, so it does not transfer any moment loading.  
+
+The choice of two BEAM189 element divisions per frame member in ANSYS is made following a mesh sensitivity study, which showed small differences between force and moment outputs for 1 and 2 element divisions, but negligible difference between results for 2 and 10 element divisions per frame member.  
+
+# 3.2. Verification results  
+
+Relative differences between the Bladed and ANSYS computed section forces, moments, displacements and rotations are shown in Table 2 for load case 1A, Table 3 for load case 1B and Table 4 for load case 2A. The difference in the peak section force values is less than $0.3\%$ for all three load cases and differences in the peak section moments between $0.0\%$ and $1.4\%$ . Where more significant percentage errors are shown for individual force and moment components, this is considered acceptable where the absolute values of these components are orders of magnitude smaller than the peak values.  
+
+Differences in the Bladed and ANSYS peak displacements and rotations for load case 1A agree to within $0.7\%$ and $0.6\%$ , respectively. For vertices where the absolute value of displacement and rotation are smaller than the peak, but not insignificant, differences are up to $5.6\%$ and $7.0\%$ , respectively. For load case 1B, differences at the peak displacement and rotations are $4.5\%$ and $2.9\%$ , respectively. Finally for load case 2A, the difference for the peak displacement is $0.7\%$ and $1.7\%$ for peak rotation.  
+
+Some discrepancies in predictions of nodal displacements and rotation are expected due to differences in the formulations for geometric non-linearity for the two software. Most importantly for the present study, however, the new method for calculation of section forces and moments in Bladed is shown to give very good agreement with ANSYS.  
+
+Table 2. Main results for load case 1A given in terms of displacements/rotations and section forces/moments at four different sections, i.e., S1 at the support, S2 at the right vertex, S3 at the left vertex, and S4 at the top. All results are given in the global frame.   
+
+
+<html><body><table><tr><td></td><td></td><td colspan="3">Displacements s(m)</td><td colspan="3">Rotations (deg)</td><td colspan="3">Force components (kN)</td><td colspan="3">Moment components (kNm)</td></tr><tr><td></td><td>Tool</td><td>Ux</td><td>Uy</td><td>Uz</td><td>x</td><td>y</td><td>z</td><td>Fx</td><td>Fy</td><td>Fz</td><td>Mx</td><td>My</td><td>Mz</td></tr><tr><td rowspan="3">S1</td><td>ANSYS</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>-55,261</td><td>-14,272</td><td>-498.2</td><td>-11,461</td><td>10,495</td><td>762,603</td></tr><tr><td>Bladed</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>-55,247</td><td>-14,264</td><td>-497.8</td><td>-11,471</td><td>10,573</td><td>766,338</td></tr><tr><td>Rel. diff.</td><td></td><td>-</td><td></td><td></td><td>-</td><td></td><td>-0.0%</td><td>-0.1%</td><td>-0.1%</td><td>+0.1%</td><td>+0.7%</td><td>+0.5%</td></tr><tr><td rowspan="3">S2</td><td>ANSYS</td><td>2.52</td><td>-2.72</td><td>0.09</td><td>0.13</td><td>-0.10</td><td>-5.50</td><td>-5,261</td><td>-9,272</td><td>1.7</td><td>894</td><td>-912</td><td>-281,053</td></tr><tr><td>Bladed</td><td>2.66</td><td>-2.60</td><td>0.09</td><td>0.12</td><td>-0.11</td><td>-5.53</td><td>-5,247</td><td>-9,264</td><td>2.2</td><td>940</td><td>-932</td><td>-281,404</td></tr><tr><td>Rel.diff.</td><td>+5.6%</td><td>+4.2%</td><td>+0.8%</td><td>+7%</td><td>-7%</td><td>-0.6%</td><td>-0.3%</td><td>-0.1%</td><td>+32%</td><td>+5%</td><td>+2%</td><td>+0.1%</td></tr><tr><td rowspan="3">S3</td><td>ANSYS</td><td>-2.00</td><td>-1.86</td><td>0.00</td><td>0.06</td><td>-0.05</td><td>-5.51</td><td>-5,261</td><td>-9,273</td><td>1.7</td><td>1,835</td><td>-1,395</td><td>-11,852</td></tr><tr><td>Bladed</td><td>-1.93</td><td>-1.93</td><td>0.00</td><td>0.06</td><td>-0.06</td><td>-5.51</td><td>-5,247</td><td>-9,264</td><td>2.2</td><td>1,860</td><td>-1,390</td><td>-12,612</td></tr><tr><td>Rel.diff.</td><td>-3.6%</td><td>+4.1%</td><td></td><td>+0.4%</td><td>+21%</td><td>+0.0%</td><td>-0.3%</td><td>-0.1%</td><td>+25%</td><td>+1.3%</td><td>-0.4%</td><td>+6.4%</td></tr><tr><td rowspan="3">S4</td><td>ANSYS</td><td>0.50</td><td>-4.60</td><td>0.07</td><td>0.10</td><td>-0.07</td><td>-3.83</td><td>-5,261</td><td>-9,272</td><td>1.8</td><td>1,126</td><td>-962</td><td>146,111</td></tr><tr><td>Bladed</td><td>0.69</td><td>-4.57</td><td>0.08</td><td>0.09</td><td>-0.08</td><td>-3.85</td><td>-5,247</td><td>-9,264</td><td>2.2</td><td>1,080</td><td>-910</td><td>146,020</td></tr><tr><td>Rel.diff.</td><td>+39%</td><td>-0.7%</td><td>+15%</td><td>-3.3%</td><td>+17%</td><td>+0.3%</td><td>-0.3%</td><td>-0.1%</td><td>+22%</td><td>-4.1%</td><td>-5.4%</td><td>-0.1%</td></tr></table></body></html>  
+
+Table 3. Main results for load case 1B given in terms of displacements/rotations and section forces/moments at four different sections, i.e., S1 at the support, S2 at the right vertex, S3 at the left vertex, and S4 at the top. All results are given in the global frame.   
+
+
+<html><body><table><tr><td></td><td></td><td colspan="3">Displacements (m)</td><td colspan="3">Rotations (deg)</td><td colspan="3">Force components (kN)</td><td colspan="3">Moment components (kNm)</td></tr><tr><td></td><td>Tool</td><td>Ux</td><td>Uy</td><td>Uz</td><td>x</td><td>y</td><td>z</td><td>Fx</td><td>Fy</td><td>Fz</td><td>Mx</td><td>My</td><td>Mz</td></tr><tr><td></td><td>ANSYS</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>-68,766</td><td>-4,766</td><td>-20</td><td>-585</td><td>378</td><td>922,425</td></tr><tr><td>S1</td><td>Bladed Rel. diff.</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>-68,945</td><td>-4,827</td><td>-19</td><td>-587</td><td>384</td><td>926,881</td></tr><tr><td></td><td></td><td></td><td></td><td>-</td><td></td><td></td><td>-</td><td>+0.3%</td><td>+1.3%</td><td>-3.6%</td><td>+0.4%</td><td>+1.4%</td><td>+0.5%</td></tr><tr><td>S2</td><td>ANSYS</td><td>2.20</td><td>-2.32</td><td>0.004</td><td>0.009</td><td>-0.003</td><td>-1.08</td><td>68,776</td><td>4,769</td><td>20</td><td>61</td><td>-25</td><td>827,599</td></tr><tr><td></td><td>Bladed</td><td>2.29</td><td>-2.23</td><td>0.003</td><td>0.006</td><td>-0.002</td><td>-1.07</td><td>68,945</td><td>4,827</td><td>19</td><td>77</td><td>10</td><td>831,745</td></tr><tr><td></td><td>Rel.diff.</td><td>+4.5%</td><td>+3.9%</td><td>+22%</td><td>+29%</td><td>+16%</td><td>+0.7%</td><td>+0.2%</td><td>+1.2%</td><td>-3.6%</td><td>+27%</td><td>-140%</td><td>+0.5%</td></tr><tr><td></td><td>ANSYS</td><td>0.92</td><td>0.94</td><td>-0.005</td><td>0.009</td><td>0.004</td><td>-1.10</td><td>31,238</td><td>5,233</td><td>-20</td><td>33</td><td>13</td><td>-580,803</td></tr><tr><td>S3</td><td>Bladed</td><td>0.94</td><td>0.92</td><td>-0.004</td><td>0.007</td><td>0.004</td><td>-1.08</td><td>31,055</td><td>5,173</td><td>-19</td><td>55</td><td>-10</td><td>-577,186</td></tr><tr><td></td><td>Rel.diff.</td><td>+1.8%</td><td>-2.1%</td><td>-20%</td><td>-26%</td><td>-14%</td><td>-2.0%</td><td>-0.6%</td><td>-1.1%</td><td>-3.6%</td><td>+70%</td><td>-179%</td><td>-0.6%</td></tr><tr><td></td><td>ANSYS</td><td>3.12</td><td>-1.39</td><td>-0.001</td><td>0.012</td><td>0.001</td><td>-2.00</td><td>73,666</td><td>5,279</td><td>23</td><td>12</td><td>-150</td><td>698,553</td></tr><tr><td>S4</td><td>Bladed</td><td>3.24</td><td>-1.33</td><td>-0.001</td><td>0.010</td><td>0.001</td><td>-2.06</td><td>73,481</td><td>5,216</td><td>23</td><td>10</td><td>-110</td><td>703,899</td></tr><tr><td></td><td>Rel.diff.+3.7%</td><td></td><td>-4.4%</td><td>-30%</td><td>-22%</td><td>-20%</td><td>+2.9%</td><td>-0.3%</td><td>-1.2%</td><td>+3.1%</td><td>-17%</td><td>-26%</td><td>+0.8%</td></tr></table></body></html>  
+
+Table 4. Main results for load case 2A given in terms of displacements/rotations and section forces/moments at five different sections, i.e., S1 at the support, S2 at the right vertex, S3 at the left vertex, S4 at the top and S5 in the bar element at the right vertex. All results are given in the global frame.   
+
+
+<html><body><table><tr><td></td><td></td><td colspan="3">Displacements (m)</td><td colspan="3">Rotations (deg)</td><td colspan="3">Force components (kN)</td><td colspan="3">Moment components (kNm)</td></tr><tr><td></td><td>Tool</td><td>Ux</td><td>Uy</td><td>Uz</td><td>x</td><td>@y</td><td>z</td><td>Fx</td><td>Fy</td><td>Fz</td><td>Mx</td><td>My</td><td>Mz</td></tr><tr><td></td><td>ANSYS</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>42,426</td><td>38,066</td><td>12.23</td><td>635.2</td><td>-471.4</td><td>-55,708</td></tr><tr><td>S1</td><td>Bladed Rel. diff.</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>42,426 -0.0%</td><td>38,108</td><td>12.08</td><td>636.6</td><td>-471.5</td><td>-55,287 -0.8%</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>+0.1%</td><td>-1.3%</td><td>+0.2%</td><td>+0.0%</td><td></td></tr><tr><td>S2</td><td>ANSYS Bladed</td><td>-0.105 -0.105</td><td>0.036</td><td>-0.005</td><td>-0.010</td><td>0.005</td><td>-0.335</td><td>42,427</td><td>-22,279</td><td>22.63</td><td>73.5</td><td>109.7</td><td>80,646 79,775</td></tr><tr><td></td><td>Rel.diff.</td><td>-0.2%</td><td>0.036 +0.2%</td><td>-0.005 -1.0%</td><td>-0.010</td><td>0.005</td><td>-0.329</td><td>42,426 -0.0%</td><td>-22,214 -0.3%</td><td>22.43</td><td>77.3</td><td>108.3 -1.2%</td><td>-1.1%</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>-1.9%</td><td>-0.4%</td><td>+1.7%</td><td></td><td></td><td>-0.9%</td><td>+5%</td><td></td><td></td></tr><tr><td>S3</td><td>ANSYS</td><td>-0.105</td><td>-0.036</td><td>0.005</td><td>-0.010</td><td>-0.005</td><td>0.335</td><td>-42,426</td><td>38,065</td><td>12</td><td>-74</td><td>110</td><td>80,643</td></tr><tr><td></td><td>Bladed</td><td>-0.105</td><td>-0.036</td><td>0.005</td><td>-0.010</td><td>-0.005</td><td>0.329</td><td>-42,426</td><td>38,108</td><td>12</td><td>-77</td><td>108</td><td>79,775</td></tr><tr><td></td><td>Rel. diff.</td><td>-0.2%</td><td>-0.2%</td><td>+1.3%</td><td>+1.9%</td><td>+0.4%</td><td>-1.7%</td><td>+0.0%</td><td>+0.1%</td><td>-1.2%</td><td>+4.4%</td><td>-1.7%</td><td>-1.1%</td></tr><tr><td>S4</td><td>ANSYS</td><td>-0.179</td><td>0.000</td><td>0.000</td><td>-0.014</td><td>0.000</td><td>0.000</td><td>0</td><td>-22,278</td><td>-20.13</td><td>0</td><td>54.6</td><td>108,353</td></tr><tr><td></td><td>Bladed</td><td>-0.178</td><td>0.000</td><td>0.000</td><td>-0.014</td><td>0.000</td><td>0.000</td><td>0</td><td>-22,214</td><td>-19.99</td><td>0</td><td>49.8</td><td>106,802</td></tr><tr><td></td><td>Rel. diff.</td><td>-0.7%</td><td></td><td></td><td>+0.4%</td><td></td><td></td><td></td><td>-0.3%</td><td>-0.7%</td><td></td><td>-8.7%</td><td>-1.4%</td></tr><tr><td>S5</td><td>ANSYS</td><td>-0.105</td><td>0.036</td><td>-0.005</td><td>-0.010</td><td>0.005</td><td>-0.335</td><td>0</td><td>60,344</td><td>-10.23</td><td>0</td><td>0</td><td>0</td></tr><tr><td></td><td>Bladed</td><td>-0.105</td><td>0.036</td><td>-0.005</td><td>-0.010</td><td>0.005</td><td>-0.329</td><td>0</td><td>60,322</td><td>-10.35</td><td>0</td><td>0</td><td>0</td></tr><tr><td></td><td>Rel. diff.</td><td>-0.2%</td><td>+0.2%</td><td>-1.0%</td><td>-1.0%</td><td>-0.4%</td><td>+1.7%</td><td></td><td>-0.0%</td><td>+1.2%</td><td></td><td></td><td></td></tr></table></body></html>  
+
+# 4. Case study: floating wind turbine support structure  
+
+The Offshore Code Comparison Collaboration (OC4) project aimed to improve the understanding and development of offshore wind energy technologies [14]. As part of the OC4 project, one of the designs examined was the DeepCwind floating semi-submersible system. This design is analysed at full scale to support the National Renewable Energy Laboratory's (NREL) offshore 5 MW baseline wind turbine [15]. The turbine serves as a representation of a large-scale, multi-megawatt structure. The DeepCwind floating semi-submersible system was selected for this case study as it is suitable for representing a statically indeterminate floating structure with multiple members. The hydrodynamic loads calculation is based on potential theory with viscous drag modelled by Morison’s equation to model the hydrodynamic damping [1].  
+
+In this example, the floating platform model has been modified to simulate a challenging scenario for calculating section forces, requiring the use of the new equilibrium-based method for section force calculation described in this paper. The unmodified OC4 DeepCwind validation model is shown in Figure 2(a). For this case study, the traditional cantilevering tower has been replaced with a slender tower supported by three pre-tensioned cables, as shown in Figure 2(b). The modified tower has a diameter of $3\textrm{m}$ and walls that are $19\,\mathrm{mm}$ thick. The cables are modelled using a massless material and are pre-tensioned to $13\;\mathrm{MN}$ in each cable. Additional point masses are added to the tower nodes to maintain the total mass and the mass centre of the original design.  The tower height interface sections remain the same as in the original OC4 model configurations detailed in [14]. It is important to note that this configuration has not been validated through DeepCwind scale model testing or numerically in the OC4 project. Instead, it is presented in this case study as a conceptual illustration.  
+
+![](88946715b61c81b1fad287b0311e456aeb3e8dba0428bd3678048ed765221a15.jpg)  
+Figure 2: (a) OC4 DeepCwind floater design (b) modified design with cable support (c) modal reference node position at tower base.  
+
+In Bladed, the floating platform and tower are considered as a single flexible body, which can undergo rigid body motions and linear structural deformations. For the purpose of calculation of Craig Bampton mode shapes [7][8], it is necessary to specify a node to be treated as fixed in space in order to exclude rigid body modes. This so-called 'modal reference node' is chosen to be located at the tower base, as shown in Figure 2(c). If the applied loading used for calculating the section forces is not in equilibrium, a non-zero reaction force will be generated at the modal reference node. This, in turn, leads to a discontinuity in the member forces at the sections on either side of the node. The purpose of this example is to demonstrate this discontinuity and how the new equilibrium-based method can serve as a better basis for calculating section forces compared to the conventional method.  
+
+A simulation was carried out with the turbine idling at a wind speed of $10\;\mathrm{m/s}$ and no waves. The section force output location and coordinate system at the modal reference node are shown in Figure 2(c). Table 5 presents the bending moment about z on the element ends on either side of this node. It is observed that the conventional method shows a $9.7\%$ discontinuity in bending moment in the structure, whereas the equilibrium-based method shows no discontinuity. The latter is the expected result as no loads are applied at the modal reference node.  
+
+Table 5. Bending moment just above and below the modal reference node (MRN) of the floater.   
+
+
+<html><body><table><tr><td></td><td colspan="3">Bending moment</td></tr><tr><td>Method</td><td>Location BelowMRN</td><td>about z (MNm) 0.964</td><td>Rel. diff. (%)</td></tr><tr><td rowspan="2">Conventional</td><td>AboveMRN</td><td>1.057</td><td rowspan="2">9.7%</td></tr><tr><td></td><td></td></tr><tr><td rowspan="2">Equilibrium</td><td>BelowMRN</td><td>1.052</td><td rowspan="2">0.0%</td></tr><tr><td>Above MRN</td><td>1.052</td></tr></table></body></html>  
+
+# 5. Conclusion  
+
+A new equilibrium-based approach for calculating internal section forces in flexible bodies has been developed and presented.  
+
+The method was used in the wind turbine design tool Bladed as an integrated part of the multibody system modelling the complete wind turbine structure. A separate method for calculating the section forces is mandatory for this application, where the deformation of the flexible bodies is described by a truncated modal expansion with a relatively low number of mode shapes (or none). The method is particularly relevant for statically indeterminate structures with significant geometrical non-linear effects, such as floating wind turbine support structure with pre-tensioned cables, where conventional methods can yield inaccuracies of the calculated section forces. In particular, this improves the confidence of practical wind turbine structural analyses in case a substructure is subject to large point loads, which may originate from pre-tensioned cables.  
+
+The new method has been verified by comparison against ANSYS, and the significance for section force calculation demonstrated through a case study based on a modified version of the DeepCwind floating platform.  
+
+# References  
+
+[1] DNV AS 2023 Bladed 4.15 theory documentaion (online) https://dnvgldocs.azureedge.net/Bladed%204.15%20Theory%20Manual/index.html   
+[2] Pellegrino P and Calladine C R 1986 Matrix analysis of statically and kinematically indeterminate frameworks Int. J. Solids Structures 22 pp 409-28   
+[3] Pellegrino P 1992 Structural computations with the singular value decomposition of the equilibrium matrix Int. J. Solids Structures 30 pp 3025-35   
+[4] Przemieniecki J S 1968 Theory of Matrix Structural Analysis (New York: McGraw-Hill)   
+[5] Cook R D, Malkus D S and Plesha M E Conceps and Application of Finite Element Analysis, $3^{r d}$ edition (New York:  John Wiley & Sons)   
+[6] Shabana A A 2014 Dynamics of Multibody Systems $\mathcal{I}^{t h}$ edition (Cambridge: Cambridge University Press)   
+[7] Craig Jr R R and Bampton M C C 1968 Coupling of substructures for dynamic analysis AIAA Journal 6 pp 1313-19   
+[8] Craig Jr R R 2000 Coupling of substructures for dynamic analysis: an overview AIAA/ASME/ ASCE/AHS/ASC 41st structures, structural dynamics and materials conference (Reston, VA: American Institute of Aeronautics and Astronautics)   
+[9]  Schwertassek R, Wallrapp O and Shabana A A 1999 Flexible Multibody Simulation and choice of shape functions Nonlinear Dynamics 20 pp 361-80   
+[10]  Schwertassek R, Dombrowski, S V and Wallrapp O 1999 Modal representation of stress in flexible multibody simulation Nonlinear Dynamics 20 pp 381-99   
+[11] Meijaard J P 2022 An extended modelling technique with generalized strains for flexible multibody systems Multibody System Dynamics 57 pp 133-55   
+[12]  Krenk S 2009 Non-linear Modeling and Analysis of Solids and Structures (Cambridge: Cambridge University Press)   
+[13]  ANSYS, Inc. 2013 ANSYS Mechanical APDL Element Reference and Theory Reference manuals, Release 14.5.7 (Canonsburg, PA: SAS IP)   
+[14]  Robertson A, et al 2014 Definition of the Semisubmersible Floating System for Phase II of OC4, NREL/TP-5000-60601 (Golden, CO: National Renewable Energy Laboratory)   
+[15]  Jonkman J, Butterfield S, Musial W and Scott G 2009 Definition of a 5-MW Reference Wind Turbine for Offshore System Development, NREL/TP-500-38060 (Golden, CO: National Renewable Energy Laboratory)  

工作OKRs/25.6-8 OKR.canvas

@@ -5,7 +5,7 @@
 		{"id":"d2c5e076ba6cf7d7","type":"text","text":"# 推进计划\n未来四周计划推进的重要事情\n\n文献调研启动\n\n建模重新推导\n\n\n","x":-600,"y":-306,"width":456,"height":347},
 		{"id":"82708a439812fdc7","type":"text","text":"# 7月已完成\n\n","x":-220,"y":134,"width":440,"height":560},
 		{"id":"505acb3e6b119076","type":"text","text":"# 6月已完成\n\n\nP1 结果对比\n- Herowind 带3.5气动与fast3.5对比 相同\n- Herowind 带4.0气动与fast4.0对比 相同\n- Herowind 带hrl气动与fast对比 需气动支持15MW\n\nP1 Bladed交流问题汇总","x":-700,"y":134,"width":440,"height":560},
-		{"id":"c18d25521d773705","type":"text","text":"# 计划\n这周要做的3~5件重要的事情，这些事情能有效推进实现OKR。\n\nP1 必须做。P2 应该做\n\n\nP1 柔性部件 叶片、塔架主动力惯性力算法 主线\n- 变形体动力学 简略看看ing\n- 柔性梁弯曲变形振动学习，主线 \n\t- 广义质量 刚度矩阵及含义\n- 如何静力学求解\n\t- 基于本构方程 读孟的论文\n\t- normal mode shape 能否使用？\n\t- F = kx 外载与弹性势能相等\n\nP1 模型线性化原理 \n- Bladed 线性化理论手册已读\n- 梳理Bladed线性化方法框架\n-\nP1 结果对比\n- Bladed与FAST之间的对比 去掉角度 不能直接比较\n- 坐标系转换\n\nP2 如何优雅的存储、输出结果。\nP1 推进气动、控制、多体、水动 耦合计算\nP2 yaw 自由度再bug确认 已知原理了\n","x":-594,"y":-693,"width":450,"height":347},
+		{"id":"c18d25521d773705","type":"text","text":"# 计划\n这周要做的3~5件重要的事情，这些事情能有效推进实现OKR。\n\nP1 必须做。P2 应该做\n\n\nP1 柔性部件 叶片、塔架主动力惯性力算法 主线\n- 变形体动力学 简略看看ing\n- 柔性梁弯曲变形振动学习，主线 \n\t- 广义质量 刚度矩阵及含义\n- 如何静力学求解 \n\t- 基于本构方程 读孟的论文\n\t- normal mode shape 能否使用？\n\t- F = kx 外载与弹性势能相等\n\t\n- 梳理bladed动力学框架 this week\n\t- 子结构文献阅读\n\nP1 模型线性化原理 this week\n- Bladed 线性化理论手册 仔细阅读\n- 梳理Bladed线性化方法框架\n-\nP1 结果对比\n- Bladed与FAST之间的对比 去掉角度 不能直接比较\n- 坐标系转换 this week\n\nP2 如何优雅的存储、输出结果。\nP2 yaw 自由度再bug确认 已知原理了\n","x":-594,"y":-693,"width":450,"height":347},
 		{"id":"30cb7486dc4e224c","type":"text","text":"# 8月已完成\n\n","x":260,"y":134,"width":440,"height":560}
 	],
 	"edges":[]

线性化求解器/Bladed线性化框架.canvas

@@ -0,0 +1,105 @@
+{
+	"nodes":[
+		{"id":"2144fa1ef85968d0","type":"text","text":"气弹模型","x":-140,"y":-400,"width":250,"height":60},
+		{"id":"bde0f2550987a8cc","type":"text","text":"线性系统","x":-140,"y":-240,"width":250,"height":60},
+		{"id":"d7b8720eb69baf7d","type":"text","text":"线性系统方程组以状态空间方程形式表示","x":-265,"y":-40,"width":250,"height":60},
+		{"id":"62dcdf23e738af0d","type":"text","text":"$$\n\\underline{{\\mathbf{x}}}=\\mathbf{x}-\\mathbf{x_{0}},\\quad\\underline{{\\mathbf{y}}}=\\mathbf{y}-\\mathbf{y_{0}},\\quad\\mathrm{~and~}\\quad\\underline{{\\mathbf{u}}}=\\mathbf{u}-\\mathbf{u_{0}}\n$$","x":85,"y":80,"width":435,"height":60},
+		{"id":"5797579dbe07ce53","type":"text","text":" $\\mathbf{x}$ 是状态向量，代表系统状态；$\\mathbf{u}$ 是系统输入向量；$\\mathbf{y}$ 是系统输出向量。","x":85,"y":180,"width":435,"height":60},
+		{"id":"58c8fd4b2b79d5c4","type":"text","text":"归一化向量 $\\underline{{\\mathbf{x}}}、\\underline{{\\mathbf{y}}}$ 和 $\\underline{{\\mathbf{u}}}$ 代表偏离平衡态的量。","x":85,"y":280,"width":435,"height":60},
+		{"id":"75f03285489bc502","type":"text","text":"矩阵 $\\mathbf{A}$、$\\mathbf{B}$、$\\mathbf{C}$ 和 $\\mathbf{D}$ 代表这些向量之间的线性化关系","x":85,"y":380,"width":435,"height":60},
+		{"id":"172345b7932e4879","type":"text","text":"$$\n\\begin{array}{r}{\\underline{{\\dot{\\mathbf{x}}}}=\\mathbf{A}\\underline{{\\mathbf{x}}}+\\mathbf{B}\\underline{{\\mathbf{u}}}}\\\\ {\\underline{{\\mathbf{y}}}=\\mathbf{C}\\underline{{\\mathbf{x}}}+\\mathbf{D}\\underline{{\\mathbf{u}}}}\\end{array}\n$$","x":178,"y":-60,"width":250,"height":100},
+		{"id":"7f80af6dea4a5833","type":"text","text":"$$\n\\begin{array}{r}{\\dot{\\mathbf{x}}=f(t,\\mathbf{x},\\mathbf{u})}\\\\ {\\mathbf{y}=h(t,\\mathbf{x},\\mathbf{u}).}\\end{array}\n$$","x":178,"y":480,"width":250,"height":100},
+		{"id":"a1a2472e2f6a21f7","type":"text","text":"弹性动力学 (Elastodynamic)：这些状态代表系统的结构模态。**弹性动力学模态由二阶运动方程控制**。因此，为了在状态空间形式中表示，每个模态由两个状态表示——位移和速度。这还包括主要的风轮刚体自由度。","x":660,"y":160,"width":283,"height":256},
+		{"id":"fe565c4907cb3e77","type":"text","text":"状态分为两大类","x":840,"y":-40,"width":250,"height":60},
+		{"id":"8d024129778a2f9a","type":"text","text":"气动 (Aerodynamic)：这些主要用于模拟动态失速和动态尾流。**这些状态通常是一阶的**，因为它们与时滞有关。","x":1000,"y":160,"width":250,"height":256},
+		{"id":"cfe15ef44610d3d4","type":"text","text":"Bladed线性化分析过程","x":-264,"y":760,"width":250,"height":60,"color":"1"},
+		{"id":"f9c41a29bddeb8cc","type":"text","text":"依次取每个工作点","x":120,"y":760,"width":250,"height":60},
+		{"id":"fb99badcdb47fb3c","type":"text","text":"找到风电机组的稳态条件，这意味着风轮没有加速，模态变形使得弹性载荷平衡外部载荷。","x":480,"y":730,"width":250,"height":120},
+		{"id":"7d322544cd4c0e9e","type":"text","text":"对于每个输入或状态，Bladed然后在平衡点两侧进行一系列幅度**逐渐增加的扰动**。**人为地增加或减少状态或输入的数值，用这些修改后的数值求解系统**，并记录状态导数和输出。扰动的数量和最大扰动幅度可以由用户定义。","x":880,"y":650,"width":250,"height":280},
+		{"id":"7798ef385eff4f6c","type":"text","text":"通过对状态导数与所有扰动值及其平衡值之间的关系进行线性回归，来**推导矩阵A、B、C和D的元素**。线性回归的梯度给出元素的值。如果相关系数小于最小相关系数，则认为该关系无效，并给该元素赋予零值。","x":1260,"y":650,"width":250,"height":280},
+		{"id":"d528580e06b6b8a1","type":"text","text":"这定义了$\\mathbf{x_{0}},\\,\\mathbf{y_{0}}$和$\\mathbf{u_{0}}$的值，即所有扰动围绕的主要平衡点。","x":465,"y":960,"width":280,"height":90},
+		{"id":"bc2d79101ffc2c67","type":"text","text":"- 需移除方位角依赖性，包括风切变、偏航、风轮不平衡等。\n- 需移除无法线性化的物理效应，例如风湍流、黏滞滑动等。","x":-302,"y":180,"width":325,"height":140},
+		{"id":"e7f67865a05e0be3","type":"text","text":"多叶片坐标系转化非旋转坐标系","x":-264,"y":1140,"width":250,"height":60},
+		{"id":"d9a4691cf95d96c3","type":"text","text":"基于 (Bir, 2008) and (Hansen, 2003)","x":-264,"y":1260,"width":250,"height":60},
+		{"id":"8d7c49775076894a","type":"text","text":"旋转坐标系到非旋转坐标系变换关系","x":465,"y":1140,"width":250,"height":60},
+		{"id":"081fcfd3683a8ce4","type":"text","text":"$$\n\\begin{align}\n\\begin{bmatrix}\nq_{0} \\\\\nq_{c} \\\\\nq_{s} \\\\\n\\end{bmatrix} &=  \n\\frac{1}{3}\\begin{bmatrix}\n1 & 1 & 1 \\\\\n2\\cos\\psi_{1} & 2\\cos\\psi_{2} & 2\\cos\\psi_{3} \\\\\n2\\sin\\psi_{1} & 2\\sin\\psi_{2} & 2\\sin\\psi_{3} \\\\\n\\end{bmatrix}\n\n\\begin{bmatrix}\nq_{1} \\\\\nq_{2} \\\\\nq_{3} \\\\\n\\end{bmatrix},\n\\end{align}\n$$","x":805,"y":1108,"width":400,"height":125},
+		{"id":"975926818ccaa396","type":"text","text":"方位角表示：$\\psi_{i}=\\Omega t+\\Psi_{i}$","x":1260,"y":1141,"width":250,"height":59},
+		{"id":"c1c28df14fb3bfbb","type":"text","text":"变换矩阵求时间导数","x":1560,"y":1140,"width":250,"height":60},
+		{"id":"edd4cdf7e36d0cfc","type":"text","text":"$$\n\\begin{align*}\n\\dot{{t}}_{R \\rightarrow NR} =\n\\frac{\\Omega}{3}\n\\begin{bmatrix}\n0 & 0 & 0 \\\\\n- 2\\sin\\psi_{1} & - 2\\sin\\psi_{2} & - 2\\sin\\psi_{3} \\\\\n2\\cos\\psi_{1} & 2\\cos\\psi_{2} & 2\\cos\\psi_{3} \\\\\n\\end{bmatrix}\n\\end{align*}\n$$","x":1880,"y":1005,"width":368,"height":135},
+		{"id":"8622b2d75c26d18b","type":"text","text":"$$\n\\begin{align*}\n\\ddot{{t}}_{R \\rightarrow NR} =\n- \\frac{\\Omega^{2}}{3}\n\\begin{bmatrix}\n0 & 0 & 0 \\\\\n2\\cos\\psi_{1} & 2\\cos\\psi_{2} & 2\\cos\\psi_{3} \\\\\n2\\sin\\psi_{1} & 2\\sin\\psi_{2} & 2\\sin\\psi_{3} \\\\\n\\end{bmatrix}\n\\end{align*}\n$$","x":1880,"y":1200,"width":368,"height":137},
+		{"id":"90ebf789028488ba","type":"text","text":"方位角带入到变换矩阵中了","x":1555,"y":1269,"width":260,"height":60},
+		{"id":"2e0560ed3e0b3d1e","type":"text","text":"系统变换矩阵","x":-264,"y":1440,"width":250,"height":60},
+		{"id":"9e5fe38ea39cb6a4","type":"text","text":"$$\n{\\bf q}_{N R}={\\bf t}_{R\\rightarrow N R}{\\bf q}_{R}\n$$","x":880,"y":1260,"width":250,"height":60},
+		{"id":"b4199d2560e40041","type":"text","text":"a common transformation matrix $\\mathbf{T}$","x":880,"y":1440,"width":250,"height":60},
+		{"id":"6980197894c4c44e","type":"text","text":"$$\n\\begin{array}{r}\n{\\bf q}_{N R}={\\bf t}_{R\\rightarrow N R}{\\bf q}_{R}\\\\\n\\dot{\\bf q}_{N R}={\\bf t}_{R\\rightarrow N R}\\dot{\\bf q}_{R}+\\dot{\\bf t}_{R\\rightarrow N R}{\\bf q}_{R}\n\\end{array}\n$$","x":408,"y":1420,"width":395,"height":100},
+		{"id":"d2862b9b5b4f20bd","type":"text","text":"$$\n\\begin{align}\n\n{T} &:= \\begin{bmatrix}\n\n{t}_{R \\rightarrow NR} & 0 \\\\\n\n\\dot{{t}}_{R \\rightarrow NR} & {t}_{R \\rightarrow NR} \\\\\n\n\\end{bmatrix}\n\n\\end{align}\n$$","x":1205,"y":1420,"width":250,"height":100},
+		{"id":"3bc2446646bd3d8b","type":"text","text":"$$\n\\begin{align}\n\n\\dot{{T}}=\n\n\\begin{bmatrix}\n\n\\dot{{t}}_{R \\rightarrow NR} & 0 \\\\\n\n\\ddot{{t}}_{R \\rightarrow NR} & \\dot{{t}}_{R \\rightarrow NR} \\\\\n\n\\end{bmatrix}.\n\n\\end{align}\n$$","x":1555,"y":1420,"width":250,"height":100},
+		{"id":"eaa37292565b4f12","type":"text","text":"$$\n\\begin{align}\n\\begin{bmatrix}\n    q_{NR} \\\\\n    \\dot{q}_{NR}\n\\end{bmatrix}\n&=\n\\begin{bmatrix}\n    t_{R \\rightarrow NR} & 0 \\\\\n    \\dot{t}_{R \\rightarrow NR} & t_{R \\rightarrow NR}\n\\end{bmatrix}\n\\begin{bmatrix}\n    q_{R} \\\\\n    \\dot{q}_{R}\n\\end{bmatrix}.\n\\end{align}\n$$","x":1880,"y":1420,"width":286,"height":100},
+		{"id":"9165677a568e0e8c","type":"text","text":"计算耦合模态","x":-264,"y":1800,"width":250,"height":60},
+		{"id":"1d1483f33284f3fb","type":"text","text":"一个包含位移和速度状态的整个状态列表的变换矩阵","x":80,"y":1420,"width":250,"height":100},
+		{"id":"9152c2d34d6164f1","type":"text","text":"坎贝尔图和叶片稳定性分析都是对**指定工作点的矩阵A**的分析","x":80,"y":1780,"width":250,"height":100},
+		{"id":"67680eb011c14d65","type":"text","text":"得到基于非旋转坐标系的模态","x":80,"y":1140,"width":250,"height":60},
+		{"id":"fad223ea1c18fca6","type":"text","text":"三叶片模型消除方位角的影响","x":80,"y":1260,"width":250,"height":60},
+		{"id":"1f219abcb5ce249f","type":"text","text":"**每个耦合模态对应一个特征值及其特征向量**","x":465,"y":1800,"width":250,"height":60},
+		{"id":"6115769fc1310020","type":"text","text":"**给定矩阵 A 的一个（复数）特征值 λ**","x":880,"y":1800,"width":250,"height":60},
+		{"id":"2203401c5015c9b8","type":"text","text":"Bladed 会根据Argand Diagram得到无阻尼频率、阻尼频率、阻尼比","x":1205,"y":1780,"width":250,"height":100},
+		{"id":"57346d11b4b1a3d0","type":"text","text":"考虑旋转坐标系下的线性模型方程，定义$\\begin{align}{x}_{R} :=\\begin{bmatrix}{q}_{R} \\\\ \\dot{{q}}_{R}\\end{bmatrix}\\end{align}$","x":1205,"y":1600,"width":305,"height":100},
+		{"id":"a680b9a9a3275ad0","type":"text","text":"矩阵A从旋转坐标系转换到非旋转坐标系","x":880,"y":1620,"width":250,"height":60},
+		{"id":"c644b36ef5f053f2","type":"text","text":"$$\n\\begin{array}{r}{\\dot{\\mathbf{x}}_{R}=\\mathbf{A}_{R}\\mathbf{x}_{R}+\\mathbf{B}_{R}\\mathbf{u}}\\\\ {\\mathbf{y}=\\mathbf{C}_{R}\\mathbf{x}_{R}+\\mathbf{D}_{R}\\mathbf{u}}\\end{array}\n$$","x":1555,"y":1595,"width":250,"height":110},
+		{"id":"32a0144ad0800646","type":"text","text":"状态向量从旋转坐标系变换到非旋转坐标系，$\\mathbf{x}_{N R}=\\mathbf{T}\\mathbf{x}_{R}$","x":1880,"y":1593,"width":250,"height":115},
+		{"id":"00f636edd7d9bb41","type":"text","text":"$$\n\\dot{{\\bf x}}_{N R}={\\bf T}\\dot{\\bf x}_{R}+\\dot{\\bf T}{\\bf x}_{R}.\n$$","x":2200,"y":1607,"width":250,"height":88},
+		{"id":"cdda60f3f9d6cdf7","type":"text","text":"$$\n\\begin{array}{r l}&{\\dot{\\mathbf{x}}_{N R}=\\mathbf{T}\\left(\\mathbf{A}_{R}\\mathbf{x}_{R}+\\mathbf{B}_{R}\\mathbf{u}\\right)+\\dot{\\mathbf{T}}\\mathbf{x}_{R}}\\\\ &{\\qquad=\\Big(\\mathbf{T}\\mathbf{A}_{R}+\\dot{\\mathbf{T}}\\Big)\\mathbf{x}_{R}+\\mathbf{T}\\mathbf{B}_{R}\\mathbf{u}}\\\\ &{\\qquad=\\Big(\\mathbf{T}\\mathbf{A}_{R}+\\dot{\\mathbf{T}}\\Big)\\mathbf{T}^{-1}\\mathbf{x}_{N R}+\\mathbf{T}\\mathbf{B}_{R}\\mathbf{u}}\\end{array}\n$$","x":2520,"y":1574,"width":334,"height":154},
+		{"id":"74b2f2439c0340f6","type":"text","text":"$$\n\\mathbf{C}_{N R}:=\\mathbf{C}_{R}\\mathbf{T}^{-1},\\quad\\mathrm{~and~}\\mathbf{D}_{N R}:=\\mathbf{D}_{R}\n$$","x":2940,"y":1678,"width":317,"height":82},
+		{"id":"e707ffb73e63911e","type":"text","text":"$$\n\\begin{array}{r l}&{\\mathbf{A}_{N R}=\\Big(\\mathbf{TA}_{R}+\\dot{\\mathbf{T}}\\Big)\\mathbf{T}^{-1}}\\\\ &{\\mathbf{B}_{N R}=\\mathbf{TB}_{R}}\\end{array}\n$$","x":2940,"y":1525,"width":250,"height":126},
+		{"id":"e1785c07b66dee47","type":"text","text":"命名，连接","x":-264,"y":2000,"width":250,"height":60},
+		{"id":"5d7b2c4818b9df8a","type":"text","text":"耦合模态如何由未耦合模态得到","x":1560,"y":1840,"width":250,"height":60,"color":"2"},
+		{"id":"849dc73be675038b","type":"text","text":"跟耦合、未耦合模态有什么关系","x":1560,"y":1760,"width":250,"height":60,"color":"2"}
+	],
+	"edges":[
+		{"id":"2970ae52ae4d1e34","fromNode":"2144fa1ef85968d0","fromSide":"bottom","toNode":"bde0f2550987a8cc","toSide":"top","label":"简化为工作点上的"},
+		{"id":"4585c4d8d627884c","fromNode":"d7b8720eb69baf7d","fromSide":"right","toNode":"172345b7932e4879","toSide":"left"},
+		{"id":"546a792fd04c1ab2","fromNode":"172345b7932e4879","fromSide":"bottom","toNode":"62dcdf23e738af0d","toSide":"top"},
+		{"id":"292e432dc4689b7b","fromNode":"62dcdf23e738af0d","fromSide":"bottom","toNode":"5797579dbe07ce53","toSide":"top"},
+		{"id":"da921e225e88c08e","fromNode":"5797579dbe07ce53","fromSide":"bottom","toNode":"58c8fd4b2b79d5c4","toSide":"top"},
+		{"id":"7b5c15166e976d55","fromNode":"58c8fd4b2b79d5c4","fromSide":"bottom","toNode":"75f03285489bc502","toSide":"top"},
+		{"id":"db2e5474468e0ab2","fromNode":"75f03285489bc502","fromSide":"bottom","toNode":"7f80af6dea4a5833","toSide":"top"},
+		{"id":"df8888162abac473","fromNode":"5797579dbe07ce53","fromSide":"right","toNode":"fe565c4907cb3e77","toSide":"left"},
+		{"id":"11267eea8de57803","fromNode":"fe565c4907cb3e77","fromSide":"bottom","toNode":"a1a2472e2f6a21f7","toSide":"top"},
+		{"id":"5e387512bde2e25b","fromNode":"fe565c4907cb3e77","fromSide":"bottom","toNode":"8d024129778a2f9a","toSide":"top"},
+		{"id":"f0fb1fb957e87422","fromNode":"cfe15ef44610d3d4","fromSide":"right","toNode":"f9c41a29bddeb8cc","toSide":"left"},
+		{"id":"71054be1129f3a45","fromNode":"f9c41a29bddeb8cc","fromSide":"right","toNode":"fb99badcdb47fb3c","toSide":"left"},
+		{"id":"02d7030b47983abe","fromNode":"fb99badcdb47fb3c","fromSide":"right","toNode":"7d322544cd4c0e9e","toSide":"left"},
+		{"id":"3d909d90f003fac6","fromNode":"7d322544cd4c0e9e","fromSide":"right","toNode":"7798ef385eff4f6c","toSide":"left"},
+		{"id":"8465ae7e85ed7cf8","fromNode":"fb99badcdb47fb3c","fromSide":"bottom","toNode":"d528580e06b6b8a1","toSide":"top"},
+		{"id":"3a2a213603dfa703","fromNode":"e7f67865a05e0be3","fromSide":"right","toNode":"67680eb011c14d65","toSide":"left"},
+		{"id":"16f797ca929927b0","fromNode":"67680eb011c14d65","fromSide":"bottom","toNode":"fad223ea1c18fca6","toSide":"top"},
+		{"id":"6b48076aa19422d2","fromNode":"67680eb011c14d65","fromSide":"right","toNode":"8d7c49775076894a","toSide":"left"},
+		{"id":"49a8654ac338efb3","fromNode":"8d7c49775076894a","fromSide":"right","toNode":"081fcfd3683a8ce4","toSide":"left"},
+		{"id":"74b3e06a95aa1d96","fromNode":"081fcfd3683a8ce4","fromSide":"right","toNode":"975926818ccaa396","toSide":"left"},
+		{"id":"02422eb8b610f951","fromNode":"c1c28df14fb3bfbb","fromSide":"right","toNode":"edd4cdf7e36d0cfc","toSide":"left"},
+		{"id":"cfe114ca82359892","fromNode":"c1c28df14fb3bfbb","fromSide":"right","toNode":"8622b2d75c26d18b","toSide":"left"},
+		{"id":"bd6a88a9a133f105","fromNode":"975926818ccaa396","fromSide":"right","toNode":"c1c28df14fb3bfbb","toSide":"left"},
+		{"id":"afc2eb67d7212673","fromNode":"c1c28df14fb3bfbb","fromSide":"bottom","toNode":"90ebf789028488ba","toSide":"top"},
+		{"id":"458c1607eb5a0e69","fromNode":"e7f67865a05e0be3","fromSide":"bottom","toNode":"d9a4691cf95d96c3","toSide":"top"},
+		{"id":"61b29bd73a9904d8","fromNode":"081fcfd3683a8ce4","fromSide":"bottom","toNode":"9e5fe38ea39cb6a4","toSide":"top"},
+		{"id":"73c967849a9e8e24","fromNode":"2e0560ed3e0b3d1e","fromSide":"right","toNode":"1d1483f33284f3fb","toSide":"left"},
+		{"id":"48de7074bd1d5808","fromNode":"1d1483f33284f3fb","fromSide":"right","toNode":"6980197894c4c44e","toSide":"left"},
+		{"id":"7c1ab81ae4bebb00","fromNode":"6980197894c4c44e","fromSide":"right","toNode":"b4199d2560e40041","toSide":"left"},
+		{"id":"9e0b611f2929c298","fromNode":"b4199d2560e40041","fromSide":"right","toNode":"d2862b9b5b4f20bd","toSide":"left"},
+		{"id":"274e7613f50a5db4","fromNode":"d2862b9b5b4f20bd","fromSide":"right","toNode":"3bc2446646bd3d8b","toSide":"left"},
+		{"id":"c75f1cdf6687afe1","fromNode":"3bc2446646bd3d8b","fromSide":"right","toNode":"eaa37292565b4f12","toSide":"left"},
+		{"id":"964767dd1d9fd76a","fromNode":"9165677a568e0e8c","fromSide":"right","toNode":"9152c2d34d6164f1","toSide":"left"},
+		{"id":"f427e9f4cb16d96e","fromNode":"9152c2d34d6164f1","fromSide":"right","toNode":"1f219abcb5ce249f","toSide":"left"},
+		{"id":"442de94eda413d36","fromNode":"1f219abcb5ce249f","fromSide":"right","toNode":"6115769fc1310020","toSide":"left"},
+		{"id":"4b63e4b900947c91","fromNode":"6115769fc1310020","fromSide":"right","toNode":"2203401c5015c9b8","toSide":"left"},
+		{"id":"2e9790fd71f8c4f8","fromNode":"2203401c5015c9b8","fromSide":"right","toNode":"849dc73be675038b","toSide":"left"},
+		{"id":"8e38d19d88116e12","fromNode":"2203401c5015c9b8","fromSide":"right","toNode":"5d7b2c4818b9df8a","toSide":"left"},
+		{"id":"3db4c58a47ae919b","fromNode":"a680b9a9a3275ad0","fromSide":"right","toNode":"57346d11b4b1a3d0","toSide":"left"},
+		{"id":"9af0ed1e8c918cb7","fromNode":"57346d11b4b1a3d0","fromSide":"right","toNode":"c644b36ef5f053f2","toSide":"left"},
+		{"id":"c176c5f030bb064e","fromNode":"32a0144ad0800646","fromSide":"right","toNode":"00f636edd7d9bb41","toSide":"left"},
+		{"id":"3d84160fe83c5ba8","fromNode":"00f636edd7d9bb41","fromSide":"right","toNode":"cdda60f3f9d6cdf7","toSide":"left"},
+		{"id":"9492ed741dc27f48","fromNode":"cdda60f3f9d6cdf7","fromSide":"right","toNode":"e707ffb73e63911e","toSide":"left"},
+		{"id":"552105f118c3ed99","fromNode":"cdda60f3f9d6cdf7","fromSide":"right","toNode":"74b2f2439c0340f6","toSide":"left"},
+		{"id":"326b7013547857dc","fromNode":"6980197894c4c44e","fromSide":"right","toNode":"a680b9a9a3275ad0","toSide":"left"}
+	]
+}

线性化求解器/参考文献/Bir-Multi-blade coordinate transformation and/auto/Bir-Multi-blade coordinate transformation and.md

@@ -0,0 +1,342 @@
+# Multi-Blade Coordinate Transformation and Its Application to Wind Turbine Analysis  
+
+Gunjit Bir1 National Renewable Energy Laboratory, Golden, Colorado, 80401  
+
+The dynamics of wind turbine rotor blades are generally expressed in rotating frames attached to the individual blades. The rotor, however, responds as a whole to excitations such as aerodynamic gusts, control inputs, and tower-nacelle motion—all of which occur in a nonrotating frame. Similarly, the tower-nacelle subsystem sees the combined effect of all rotor blades, not the individual blades. Multi-blade coordinate transformation (MBC) helps integrate the dynamics of individual blades and express them in a fixed (nonrotating) frame. MBC  involves  two  steps:  transformation  of  the  rotating  degrees  of  freedom,  and transformation of the equations of motion. This paper details the MBC operation. A new MBC scheme is developed that is applicable to variable-speed turbines, which may also have dissimilar blades. The scheme also covers control, disturbance, output, and feed-forward matrices. Depending on the analysis objective, wind turbine researchers may generate system matrices either in the first-order (state-space) form or the second-order (physicaldomain) form. We develop MBC relations for both these forms. MBC is particularly essential for modal and stability analyses. Commonly, wind turbine researchers first compute the periodic state-space matrix, time-average it over the rotor rotational period, and then apply conventional eigenanalysis to compute modal and stability characteristics. Direct averaging, however, eliminates all periodic terms that contribute to system dynamics, thereby producing errors. While averaging itself is not always a bad approach, it must follow MBC. Sample results are presented to illustrate this point and also to show the application of MBC to the modal and stability analysis of a 5-MW turbine.  
+
+Nomenclature  
+
+A $=$ system dynamics matrix   
+$B$ $=$ control matrix for the first-order system   
+$B_{d}$ $=$ disturbance matrix for the first-order system   
+$b$ $=$ blade number   
+$C$ $=$ matrix consisting of damping and gyroscopic terms   
+$F$ $=$ control matrix for the second-order system   
+$F_{d}$ $=$ disturbance matrix for the second-order system   
+$K$ $=$ matrix consisting of centrifugal terms and aerodynamic and structural stiffness terms   
+$m$ $=$ the number of rotating degrees of freedom per blade   
+$M$ $=$ mass matrix   
+$n F$ $=$ nonrotating (fixed) degrees of freedom   
+$N$ $=$ number of rotor blades   
+$X$ $=$ physical vector comprising the system degrees of freedom   
+$X_{F}$ $=$ vector of degrees of freedom for the nonrotating part of the system   
+$X_{N R}$ $=$ vector of system degrees of freedom following MBC   
+$\widetilde{t_{1}},\widetilde{t_{2}},\widetilde{t_{3}}$ = 3x3 transformation matrices, which are functions of the rotor azimuth   
+$T_{I},\,T_{2},\,T_{3}$ =  block-diagonal transformation matrices, which are functions of the rotor azimuth   
+u $=$ control vector   
+$w$ $=$ disturbance vector   
+Y $=$ output vector   
+z $=$ state-space vector   
+$q_{b}$ $=$ bth blade degree of freedom in the rotating frame   
+$\boldsymbol{\mathcal{O}}$ $=$ rotor angular speed   
+$\mathcal{V}$ $=$ reference (1st) blade azimuth   
+$\Psi_{b}$ $=$ bth  blade azimuth   
+$(~)_{N R}$ $=$ a quantity referred to in the nonrotating frame  
+
+# I.  Introduction  
+
+Th e dynamics of wind turbine rotor blades are generally expressed as rotating frames attached to the individual blades. The rotor, however, responds as a whole to excitations such as aerodynamic gusts, control inputs, and tower-nacelle motion—all of which occur in a nonrotating frame.  Similarly, the tower-nacelle subsystem sees the combined effect of all rotor blades, not the individual blades. Multi-blade coordinate transformation (MBC) helps integrate the dynamics of individual blades and express it in a fixed (nonrotating) frame. MBC offers several benefits. First, it properly models the dynamic interaction between the nonrotating tower-nacelle and the spinning rotor. Second, it offers physical insight into rotor dynamics and how the rotor interacts with fixed-system entities, such as wind, controls, and tower-nacelle subsystem. Third, MBC filters out all periodic terms except those which are integral multiples of $\Omega\mathbf{N}$ , where $\Omega$ is the rotor angular speed and N is the number of rotor blades. A wind turbine system is basically a periodic system; equations governing its dynamics show periodic parameters, which arise because of the periodic interaction between the rotating subsystem (rotor) and the nonrotating entities (tower, nacelle, wind, controls, and gravity). The blade equations usually contain all harmonics. MBC shows that a rotor behaves as a filter, allowing only specific harmonics of blade motion to be felt by the fixed system. This filtering action also renders the system equations numerically well-conditioned; all nonessential periodic terms are eliminated.  
+
+MBC is widely used in the helicopter field. Miller used it to analyze the flap-motion related stability and control. Coleman and Feingold2 used it to analyze the rotor in-plane motion (lag motion); it was the first successful attempt to understand the helicopter ground resonance problem, which had been eluding the earlier researchers. However, these efforts applied the MBC only in a heuristic fashion. Hohenemser and $\mathrm{Yin}^{3}$ provided the first mathematically sound basis. Later, Johnson provided a systematic mathematical basis and thorough exposition of the MBC. Using this mathematical basis, Bir et al5 developed a numerical MBC approach that could be applied to a general helicopter system governed by arbitrary degrees of freedom; the approach was used for stability and response analysis of several helicopters.  
+
+Because MBC offers so many benefits, it is receiving attention in the wind turbine field. Bir and Butterfield included MBC in a stability analysis scheme and predicted an instability caused by coalescence between the rotor inplane and the tower motions. They also showed how MBC provides a physical insight into the rotor in-plane motion. In the turbine field, Malcolm appears to be the first to provide a mathematical form of the turbine equations following application of the MBC.  His prime motivation was to relate the inflow characteristics with the turbine response and to extract linearized models from general aeroelastic codes such as ADAMS. McCoy extended this effort to obtain linear time-invariant system equations required in the standard control design approaches. Hansen used MBC for improved modal dynamics to avoid stall-induced vibrations10 and later combined it with an eigenvalues approach to study the aeroelastic stability characteristics of a three-bladed turbine . Riziotis et al applied MBC to analyze stability of two three-bladed turbines: one stall-regulated and the other pitch-regulated. Bir and Jonkman used MBC in conjunction with FAST   to study aeroelastic characteristics of a 5-MW turbine in both land-based and offshore configurations14.  
+
+All attempts at MBC thus far, both in the helicopter and wind turbine fields, have assumed the rotor speed to be constant and the rotor blades to be similar. A modern wind turbine is rarely constant-speed. Also, turbines may not have identical blades, structurally or aerodynamically. We need an MBC approach that overcomes the two limitations.  
+
+This paper provides a new MBC scheme that is applicable to a variable-speed turbine, which may also have dissimilar blades. The scheme also covers control, disturbance, output, and feed-forward matrices, which have been ignored to date. Depending on the analysis objective, wind turbine researchers may generate system matrices either in the first-order (state-space) form or the second-order (physical-domain) form. We develop MBC relations for both these forms. In literature, MBC is also referred to as the Fourier coordinate transformation (FCT) and as the Coleman transformation.  
+
+Section II describes the MBC operation and how it relates the blade and rotor coordinates. Section III shows how these basic relations are used to transform the system equations of motion from rotating to nonrotating (fixed) coordinates, which is the objective of MBC. We present the two possible approaches to MBC transformation, the substitution method and the operational method, and explain why the substitution method is better suited to numerical analysis. To demonstrate the application of MBC, we chose a conceptual 5-MW, three-blade, upwind wind turbine15. Section IV brief describes the properties of this turbine. Though MBC has several applications, it is essential to and mostly used for modal and stability analyses. Commonly, wind turbine researchers first compute the periodic state-space matrix, time-average it over the rotor rotational period, and then apply conventional eigenanalysis to compute modal and stability characteristics. The averaging, however, eliminates all periodic terms that contribute to system dynamics and can cause erroneous results. While averaging itself is not a bad approach, it must follow MBC (Strictly, a Floquet analysis—not averaging—must follow MBC). Section V presents sample results that show how MBC is essential to correct stability analysis.  
+
+# II.  MBC: Transformation of Coordinates  
+
+Generally, the turbine equations (i.e., the coupled tower-nacelle-rotor equations) are derived using mixed degrees of freedom, some of which may be in the rotating frame and the other in the nonrotating frame. This is sometimes desirable. For example, in some simulations studies, we may be interested in studying the tower response in the ground-fixed (nonrotating) frame and the blades response in their respective rotating frames. Often, however, we are interested in understanding the coupled behavior of tower-nacelle-rotor system. In such cases, it is desirable to express the full system behavior in a fixed frame. MBC helps us achieve this through a rotating-frame to nonrotating-frame coordinate transformation.  
+
+Consider a rotor with $N$ blades that are spaced equally around the rotor azimuth. In such a case, the azimuth location of $b$ th blade is given by  
+
+$$
+\psi_{b}=\psi+(b-1)\frac{2\pi}{N}
+$$  
+
+where $\mathcal{V}$ is the azimuth of the first (reference) blade. We assume that $\psi=0$ implies the first blade is vertically up. Let $q_{b}$ be a particular rotating degree of freedom for the $b$ th blade. The MBC is a linear transformation that relates the rotating degrees of freedom to new degrees of freedom defined as follows [4]:  
+
+$$
+\begin{array}{l}{{q_{0}=\displaystyle{\frac{1}{N}\sum_{b=1}^{N}q_{b}}}}\\ {{\displaystyle q_{n c}=\displaystyle{\frac{2}{N}\sum_{b=1}^{N}q_{b}}\cos n\,{\psi_{b}}}}\\ {{\displaystyle q_{n s}=\displaystyle{\frac{2}{N}\sum_{b=1}^{N}q_{b}}\sin n\,{\psi_{b}}}}\\ {{\displaystyle q_{N_{2}}=\displaystyle{\frac{1}{N}\sum_{b=1}^{N}q_{b}(-1)^{b}}}}\end{array}
+$$  
+
+(2)  
+
+These new degrees of freedom are in the nonrotating (fixed) frame. In literature, these are called nonrotating degrees of freedom; we call these rotor coordinates because they express the cumulative behavior of all rotor blades (and not individual blades) in the fixed frame. The physical interpretation of each rotor coordinate depends on the degree of freedom it refers to. For example, if $q_{b}$ is a flap degree of freedom, then $q_{O}$ is the rotor coning, $q_{I c}$ is the rotor tippath-plane fore-aft tilt about an axis that is horizontal and normal to the rotor shaft, and $q_{I s}$ is the rotor tip-path-plane side-side tilt about an axis that is vertical and normal to the rotor shaft. If $q_{b}$ is a lag degree of freedom, then $q_{O}$ is the rotor collective lag, $q_{I c}$ is the horizontal displacement of the rotor center-of-mass in the rotor plane, and $q_{I s}$ is the vertical displacement of the rotor center-of-mass in the rotor plane [6]. Other degrees of freedom may be interpreted using equations (3) and drawing supporting sketches. The rotor modes corresponding to $q_{n c}$ and $q_{n s}\left(n\!>\!I\right)$ and $q_{N/2}$ are called reactionless modes because they do not cause any transference of moments or forces from the rotor to the hub (fixed frame). The value of n goes from $^{\,l}$ to $(N{-}I)/2$ if $N$ is odd, and from $^{\,l}$ to $(N\!-\!2)/2$ if $N$ is even. The $q_{N/2}$ mode, called the differential mode, exists only if the number of blades is even.  
+
+Most of the wind turbines in the world are three-bladed. Therefore, we will assume $N$ to be $^3$ in the rest of this paper. The formulation developed in the paper for $N{=}3$ , however, is general and can be readily extended for more blades. Setting $N{=}3$ in equations (2), we obtain  
+
+$$
+\begin{array}{l}{{q_{0}=\displaystyle{\frac{1}{3}\sum_{b=1}^{3}}q_{b}}}\\ {{\displaystyle q_{c}=q_{1c}=\frac{2}{3}\sum_{b=1}^{3}q_{b}\cos\psi_{b}}}\\ {{\displaystyle q_{s}=q_{1s}=\frac{2}{3}\sum_{b=1}^{3}q_{b}\sin\psi_{b}}}\end{array}
+$$  
+
+We will call $q_{I c}$ the cosine-cyclic mode and $q_{I s}$ the sine-cyclic mode. These two cyclic modes, together with the coning mode, $q_{O}$ , lead to coupling of the rotor with the rest of the turbine.  Equations (3) determine the rotor coordinates, given the blade coordinates. The inverse transformation, yielding the blade coordinate given the rotor coordinates, is  
+
+$$
+q_{b}=q_{0}+q_{c}\cos\psi_{b}+q_{s}\sin\psi_{b};\;\;\;b=1,2,3
+$$  
+
+# III.  MBC: Transformation of Equations of Motion  
+
+Most aeroelastic codes generate linear equations in the second-order form and, if required (for controls, for example), transform these to a first-order form. We derive MBC relations for both the second- and first-order systems.  
+
+# A. Second-Order System Equations  
+
+Second-order system equations may be written as  
+
+$$
+M\ddot{X}+C\dot{X}+K X=F u+F_{d}w
+$$  
+
+The associated output equations are  
+
+$$
+Y=C_{\nu}\dot{X}+C_{d}X+D u+D_{d}w
+$$  
+
+where $M,\,C,\,K,\,F,$ , and $F_{d}$ are respectively the mass, damping, stiffness, control, and disturbance matrices. The $u$ and $w$ are respectively the control and disturbance vectors. Note that the $\mathrm{C}$ matrix contains gyroscopic terms in addition to the structural damping terms, even though we designate it as a damping matrix. Similarly, K contains centrifugal terms in addition to the structural and aerodynamic stiffness terms. Also, each of these matrices contains direct and cross-coupling terms. For example, $M$ contains direct blade inertias as well as blade-tower coupling inertias. However, the formulation we develop does not require a user to explicitly partition any of the matrices into submatrices delineating direct and cross-coupling terms. The partitioning is implicit in the definition of the coordinates vector $X$ (also called the physical vector in literature). Assume that $X$ is expressed as  
+
+where $X_{F}$ is a $n F{\times}I$ column vector representing the $n F$ fixed-frame-referenced degrees of freedom and $q_{b}^{j}$ is the $j^{t h}$ rotating degree of freedom for the $b^{t h}$ blade. As is evident from (6), $m$ is the number of rotating degrees of freedom for each blade. Thus the length of vector $X$ is $n F{+}3m$ , the total number of degrees of freedom for the full system. Most often, aeroelastic codes assume that the physical vector has the form (6). If not, we can always perform simple rows and columns exchange to transform system equations (5) such that the physical vector $\Chi$ assumes the form (6).  
+
+There are two methods to transform equations (5) to the nonrotating frame: the operational method and the substitution method. In the operational method, we apply the following summation operators to the rotating-frame equations of motion:  
+
+$$
+\begin{array}{l}{{\displaystyle\frac{1}{3}\sum_{b=1}^{3}(r o t_{-}e q n s)}}\\ {~~}\\ {{\displaystyle\frac{2}{3}\sum_{b=1}^{3}(r o t_{-}e q n s)\cos\psi_{b}}}\\ {~~}\\ {{\displaystyle\frac{2}{3}\sum_{b=1}^{3}(r o t_{-}e q n s)\sin\psi_{b}}}\end{array}
+$$  
+
+where rot-eqns refers to the rotating-frame equations of motion. Note the similarity of summation operators in equations (3) and (7). To accomplish these operations, we first need that the rot-eqns are available in an analytical form showing all periodic terms explicitly; the periodic terms involve trigonometric functions of $\mathrm{sin}(k\,\Psi_{b})$ and $\cos(k\,\varPsi_{b})$ terms, where $\boldsymbol{\mathrm{k}}$ in general can have any integer value depending on the system. First, we must express all products of trigonometric functions as sums of trigonometric functions. Next, we must perform cumbersome operations to ensure that the value of each $n$ in these trigonometric functions satisfies the requirements of equations (3). Finally, we use transformation (3) to convert the rot-eqns to the nonrotating frame. The operational method is thus quite tedious and applicable only if we have equations of motion available in an analytical form.  
+
+All aeroelastic codes generate equations of motion in a numerical form and, therefore, are not amenable to the operational method. We must instead use the substitutional method. In this method, we substitute the rotational degrees of freedom with the rotor coordinate using equation (4). If we do so, the three rotational degrees of freedom $(q_{1}^{j},q_{2}^{j},q_{3}^{j})$ , corresponding to the jth degree of freedom for each of the three blades, may be transformed to nonrotating coordinates $(q_{0}^{j},q_{c}^{j},q_{s}^{j})$ using the following relation:  
+
+$$
+\left\{q_{1}^{j}\right\}=\widetilde{t}\left\{q_{0}^{j}\atop q_{2}^{j}\right\}=\widetilde{t}\left\{q_{c}^{j}\atop q_{s}^{j}\right\}
+$$  
+
+where  
+
+$$
+{\widetilde{t\,}}={\left[\begin{array}{l l l}{1}&{\cos\psi_{1}}&{\sin\psi_{1}}\\ {1}&{\cos\psi_{2}}&{\sin\psi_{2}}\\ {1}&{\cos\psi_{3}}&{\sin\psi_{3}}\end{array}\right]}
+$$  
+
+Using (8), the full-system rotating-frame degrees-of-freedom vector, $X$ , is expressed in terms of the nonrotatingframe degrees-of-freedom vector, $X_{N R}$ , as follows  
+
+$$
+X=T_{1}X_{N R}
+$$  
+
+where $T_{l}$ is the block diagonal matrix:  
+
+$$
+T_{1}=\left[\begin{array}{l l l l l}{I_{n F\times n F}}&&&&\\ &{\widetilde{t}}&&&\\ &&{\widetilde{t}}&&\\ &&&{\ddots}&\\ &&&&{\widetilde{t}}\end{array}\right]_{(n F+3m)},
+$$  
+
+and  
+
+$$
+\begin{array}{r}{X_{n c}=\left\lfloor\begin{array}{l}{X_{t}}\\ {q_{0}^{1}}\\ {q_{1}^{2}}\\ {q_{1}^{3}}\\ {q_{0}^{4}}\\ {q_{1}^{5}}\\ {q_{2}^{6}}\\ {q_{1}^{7}}\\ {\vdots}\\ {q_{n}^{9}}\\ {q_{n}^{10}}\end{array}\right\rfloor}\end{array}
+$$  
+
+Taking the first and second time derivatives of the two sides of equation (10), we obtain  
+
+$$
+\begin{array}{r l}&{\dot{X}=T_{1}\dot{X}_{_{N R}}+\Omega T_{2}X_{_{N R}}}\\ &{\ddot{X}_{_{N R}}=T_{1}\ddot{X}_{_{N R}}+2\Omega T_{2}\dot{X}_{_{N R}}+(\Omega^{2}T_{3}+\dot{\Omega}T_{2})X_{_{N R}}}\end{array}
+$$  
+
+where $\Omega$ is the rotor angular velocity and Ω& is the rotor angular acceleration. The $T_{2}$ and $T_{3}$ transformation matrices are given by  
+
+$$
+T_{2}=\left[\begin{array}{l l l l l}{0_{n F\times n F}}&&&&\\ &{\widetilde{t}_{2}}&&&\\ &&{\widetilde{t}_{2}}&&\\ &&&{\ddots}&\\ &&&&{\widetilde{t}_{2}}\end{array}\right]_{(n F+3)}
+$$  
+
+$$
+T_{3}=\left[\begin{array}{l l l l l}{0_{n F\times n F}}&&&&\\ &{\tilde{t}_{3}}&&&\\ &&{\tilde{t}_{3}}&&\\ &&&{\ddots}&\\ &&&&{\tilde{t}_{3}}\end{array}\right]_{(n F+3m)\times(n F+3m)}
+$$  
+
+where  
+
+$$
+\widetilde{t}_{2}=\left[\begin{array}{c c c}{0}&{-\sin\psi_{1}}&{\cos\psi_{1}}\\ {0}&{-\sin\psi_{2}}&{\cos\psi_{2}}\\ {0}&{-\sin\psi_{3}}&{\cos\psi_{3}}\end{array}\right]\quad a n d\quad\widetilde{t}_{3}=\left[\begin{array}{c c c}{0}&{-\cos\psi_{1}}&{-\sin\psi_{1}}\\ {0}&{-\cos\psi_{2}}&{-\sin\psi_{2}}\\ {0}&{-\cos\psi_{3}}&{-\sin\psi_{3}}\end{array}\right]
+$$  
+
+Substituting for $X$ and its time derivatives from equations (10) and (13) in equations (5a) and (5b), we obtain the system equations of motion in the nonrotating frame as  
+
+$$
+M_{\scriptscriptstyle{N R}}\ddot{X}_{\scriptscriptstyle{N R}}+C_{\scriptscriptstyle{N R}}\dot{X}_{\scriptscriptstyle{N R}}+K_{\scriptscriptstyle{N R}}X_{\scriptscriptstyle{N R}}=F_{\scriptscriptstyle{N R}}u_{\scriptscriptstyle{N R}}+F_{d_{\scriptscriptstyle{N R}}}w
+$$  
+
+and  
+
+$$
+Y_{_{N R}}=C_{\nu_{N R}}\dot{X}_{_{N R}}+C_{d_{N R}}X_{_{N R}}+D_{_{N R}}u_{_{N R}}+D_{_{d_{N R}}}w
+$$  
+
+The subscript NR identifies the associated quantity to be in the nonrotating frame. The various matrices appearing in equations (17) and (18) are as follows  
+
+$$
+\begin{array}{r l}&{M_{N R}=M T_{1}}\\ &{C_{N R}=2\Omega M T_{2}+C T_{1}}\\ &{K_{N R}=\Omega^{2}M T_{3}+\dot{\Omega}M T_{2}+\Omega C T_{2}+K T_{1}}\\ &{F_{N R}=F T_{1c}}\\ &{F_{d N R}=F T_{d}}\\ &{C_{\nu R}=T_{1}^{-1}C_{\nu}T_{1}}\\ &{C_{d N R}=T_{1}^{-1}(\Omega C_{\nu}T_{2}+C_{d}T_{1})}\\ &{D_{N R}=T_{10}^{-1}D T_{1c}}\\ &{D_{d N R}=T_{1}^{-1}D J_{d}}\end{array}
+$$  
+
+Note that the variable-speed operation affects only the stiffness matrix. Also, the disturbance matrix stays unchanged; this is because the disturbance $\boldsymbol{w}$ is already in the nonrotating frame. If $F c$ and mc are the number of controls in the fixed and rotating frames respectively, then the control vectors $u$ and $u_{M R}$ are related as follows:  
+
+$$
+u=T_{1c}u_{N R}
+$$  
+
+where  
+
+$$
+T_{1c}=\left[\begin{array}{l l l l l}{T_{F e x F c}}&&&&\\ &{\widetilde{t}}&&&\\ &&{\widetilde{t}}&&\\ &&&{\ddots}&\\ &&&&{\widetilde{t}}\end{array}\right]_{(F c+3m c)X(F_{4})}
+$$  
+
+Similarly,  
+
+$$
+Y=T_{1o}Y_{N R}
+$$  
+
+where  
+
+where $F o$ and mo are the number of outputs in the fixed and rotating frames respectively. In summary, following the application of MBC, the rotating-frame second-order equations (5a) and (5b) are transformed to the nonrotatingframe equations (17) and (18).  
+
+# B. First-Order System Equations  
+
+First-order system equations are generally expressed as  
+
+$$
+\dot{z}=A z+B u+B_{d}w
+$$  
+
+where $z$ is the state-space vector. The $A,B,$ , and $B_{d}$ are respectively the system dynamics, control, and disturbance matrices. The $u$ and $w$ are the control and disturbance vectors. The associated output equations are  
+
+$$
+\boldsymbol{Y}=\boldsymbol{C}\boldsymbol{z}+\boldsymbol{D}\boldsymbol{u}+\boldsymbol{D}_{d}\boldsymbol{w}
+$$  
+
+The state-space vector, $z$ , is related to the physical vector, $X$ , as follows:  
+
+$$
+z=\left\{{X\atop\dot{X}}\right\}
+$$  
+
+Using equations (26) , (8-16) and (20-23), equations (24) and (25) are transformed to the nonrotating frame as follows:  
+
+$$
+\dot{z}_{N R}=A_{N R}z_{N R}+B_{N R}u_{N R}+B_{d N R}w
+$$  
+
+and  
+
+$$
+Y_{N R}=C_{N R}z_{N R}+D_{N R}u_{N R}+D_{d_{N R}}w
+$$  
+
+The subscript NR identifies the associated quantity to be in the nonrotating frame. The various matrices appearing in equations (27) and (28) are as follows  
+
+$$
+\begin{array}{r l}&{A_{N R}=\left[\!\!\begin{array}{c c}{T_{1}^{-1}}&{0}\\ {0}&{T_{1}^{-1}\bigg]\bigg\{\boldsymbol A\left[\begin{array}{c c}{T_{1}}&{0}\\ {0T_{2}}&{T_{1}\bigg]\!-\!\left[\Omega T_{2}}&{0}\\ {0^{2}T_{3}+\hat{\Omega}T_{2}}&{2\Omega T_{2}}\end{array}\right]\!\!\right\}}\\ &{B_{N R}=\left[\!\!\begin{array}{c c}{T_{1}^{-1}}&{0}\\ {0}&{T_{1}^{-1}\bigg]B T_{1c}}\\ {C_{N R}=T_{1}^{-1}\big[C_{1}T_{1}+\Omega C_{2}T_{2}}&{C_{2}T_{1}\big];\quad C=\left[C_{1}\,}&{C_{2}\right]}\end{array}\right.}\\ &{D_{N R}=T_{10}^{-1}D T_{1C}}\\ &{D_{d N R}=T_{10}^{-1}D_{d}}\end{array}
+$$  
+
+In summary, following the application of MBC, the rotating-frame first-order equations (24) and (25) are transformed to the nonrotating-frame equations (27) and (28). Also, note that the dynamics matrix, $A$ , is influenced by the variable rotor speed.  
+
+# C. Misconceptions About MBC  
+
+Certain misconceptions regarding MBC frequently arise both in the helicopter and the wind energy fields. The main ones are:  
+
+1) The rotor must spin at a constant speed.   
+2) MBC transforms a time-variant system to a time-invariant one.   
+3) The rotor blades must have identical structural and aerodynamic properties.  
+
+None of the above statements is true. Explanatory remarks follow:  
+
+1) In the previous section, we derived MBC relations assuming that the rotor speed, $\Omega$ , varies. MBC is thus applicable to variable-speed turbines.   
+2) MBC does not eliminate all periodic terms, thereby transforming a time-variant system to a time-invariant one. MBC actually acts as a filter; it lets through terms that are integral multiples of $\Omega\mathbf{N}$ and filters out the rest the other periodic terms, provided all blades are identical and the operating conditions are timeinvariant.   
+3)  In deriving the MBC relations, we never assumed the rotor blades to be identical, structurally or aerodynamically. Any dissimilarity will be reflected in the system matrices and MBC can be applied to such matrices. For dissimilar blades, however, MBC will not be able to filter out even those periodic terms that are not integral multiples of ΩN.   
+The only restriction for MBC to be applicable is that the blades be spaced equally around the rotor azimuth.   
+Another misconception (related to the second misconception listed above) is that Floquet analysisis not required  
+
+if MBC is applied. Because MBC does not eliminate all periodic terms, a Floquet4,16 analysis is required, particularly if the blades are dissimilar.  
+
+# IV.  Brief Description of the 5-MW Wind Turbine  
+
+The conceptual wind turbine used for stability analysis in this paper consists of an upwind, three-bladed, 126-meter diameter rotor mounted on top of an 87.6-meter tower. Jonkman et al [15] selected its properties so as to best match a few existing 5-MW designs. Salient turbine properties are summarized in Table 1; other properties can be found in the cited reference. The NREL offshore 5-MW baseline wind turbine has been used to establish the reference specifications for a number of research projects supported by the U.S. Department of Energy’s Wind Energy Technologies Program.  
+
+Table 1. Summary of Baseline Wind Turbine Properties   
+
+
+<html><body><table><tr><td>Rating</td><td>5MW</td></tr><tr><td>RotorOrientation</td><td>Upwind</td></tr><tr><td>Control</td><td>VariableSpeed,CollectivePitch</td></tr><tr><td>RotorDiameter</td><td>126m</td></tr><tr><td>Hub Height</td><td>90 m</td></tr><tr><td>Cut-In, Rated, Cut-Out Wind Speed</td><td>3 m/s,11.4 m/s,25 m/s</td></tr><tr><td>Cut-ln,RatedRotorSpeed</td><td>6.9 rpm,12.1 rpm</td></tr><tr><td>RatedTipSpeed</td><td>80 m/s</td></tr><tr><td>Overhang,ShaftTilt,Precone</td><td>5m,5°,2.5°</td></tr><tr><td>BladeStructuralDampingRatio</td><td>1.0%</td></tr><tr><td>TowerStructuralDampingRatio</td><td>1.0%</td></tr></table></body></html>  
+
+# V.  Results: Application of MBC to Stability Analysis  
+
+For stability analysis, we consider two sets of operating conditions: parked conditions and normal operating conditions. In the parked condition, the generator torque is assumed to be zero, implying that the rotor is free to idle. All results are obtained using a structural damping ratio of $1\%$ for both the tower and the blades. Because the linearization scheme in FAST cannot yet handle active controls and dynamic stall, we ignore these effects during stability analyses.  
+
+# A. Parked Condition  
+
+We analyzed stability characteristics of the 5-MW turbine for a number of parked conditions, which are mandated by the International Electrotechnical Commission (IEC) design standard17,18. However, to illustrate the importance of MBC, we present results for only one parked condition covered under the IEC load case 6.2a. Due to loss of grid, the turbine is free to idle. It encounters an extreme $50\;\mathrm{m/s}$ wind with the nacelle yaw set at $140^{0}$ with respect to the wind direction. All blades are feathered at $90^{0}$ . Under these conditions, trim analysis using FAST predicts that the rotor would idle at an angular speed of $0.0076~\mathrm{Hz}$ . Using this rotor speed, we again used FAST to  
+
+generate the system dynamic matrix, A. (see equation 24). Next, we used two approaches to compute the stability   
+characteristics. In the first approach, we azimuth-averaged the system matrix $A$ and then computed its eigensolution.   
+Eigenvalues provided the modal damping and frequencies; eigenvectors helped us identify the modes. Table 2 lists   
+the modal dampings. A negative damping, identified in bold red, implies instability.  
+
+In the second approach, we first performed MBC to obtain $A_{N R}$ . Then we azimuth-averaged the matrix $A_{N R}$ and computed its eigensolution. The third column in Table 2 lists the predicted modal dampings.  
+
+Table 2. Modal Damping Ratios (nacelle yaw $\mathbf{=}\mathbf{140^{0}}$ ; all blades pitched at $90^{\boldsymbol{0}}.$ )   
+
+
+<html><body><table><tr><td rowspan="2">Mode</td><td colspan="2">Damping ratio (%)</td></tr><tr><td>WithoutMBC</td><td>WithMBC</td></tr><tr><td>1st lag</td><td>0.18</td><td>-1.48</td></tr><tr><td>1st flap</td><td>1.07</td><td>2.32</td></tr><tr><td>2nd flap</td><td>2.28</td><td>1.71</td></tr><tr><td>tower 1st side-to-side</td><td>-1.93</td><td>2.05</td></tr><tr><td>tower 1st fore-aft</td><td>0.68</td><td>2.19</td></tr><tr><td>tower 2nd side-side</td><td>3.84</td><td>2.02</td></tr><tr><td>tower 2nd fore-aft</td><td>0.96</td><td>2.87</td></tr></table></body></html>  
+
+We see that the two approaches not only yield different damping predictions, but can also lead to erroneous conclusions. If MBC is not used, instability might be interpreted as stability and vice versa. The reason is simple. In the first approach, direct averaging eliminates all periodic terms that contribute to system dynamics, therefore, erroneous results follow. In the second approach, the MBC is applied first; this appropriately transforms the periodic terms and also filters out the nonessential periodicity. Results are, therefore, improved. While averaging itself is not a bad approach, it must follow MBC (strictly, a Floquet analysis—not averaging—must follow MBC).  
+
+# B. Normal Operating Condition  
+
+As is done in a typical stability analysis, we first obtain modal frequencies versus the rotor speed when the rotor is spinning in vacuum. These frequencies and associated modes offer insight into system dynamics and help interpret the stability modes. The results are obtained using MBC followed by eigenanalysis. This plot of modal frequencies versus the rotor speed, usually called the Coleman diagram, is shown in Fig.1.  
+
+A consequence of the MBC transformation, described in Section III, is that the individual blade modes are transformed to rotor modes. Consider the first lag (edgewise) mode of each blade. Instead of seeing three first lag modes, one for each blade, we see three rotor modes: Lag 1_coll, Lag 1_prog, and Lag 1_reg. These modes are respectively the first-lag collective, progressive, and regressive modes and are physically more meaningful. Note that the tower-nacelle-drivetrain system sees the cumulative effect of all blades and not the individual blades. The rotor modes in fact represent the cumulative effect of all the blades. The collective lag mode is a rotor mode in which all the blades bend, within the rotor plane, in the same direction (either clockwise or counterclockwise) and with the same magnitude and phase. One can visualize that this rotor mode would interact with the drivetrain torsion mode. In the other two lag modes, the blades oscillate within the plane of the rotor in such a way that the effective center of mass of the rotor whirls about the rotor shaft. In the progressive mode, the center of mass whirls in the same direction as the rotor spin direction but at an angular speed higher than the rotor speed. In the regressive mode, the center of mass whirls in a direction opposite to that of the rotor spin and at an angular speed lower than the rotor speed. One can visualize that the mass imbalance associated with the rotor center-of-mass whirl can induce whirling of the rotor shaft bending mode and to some extent the tower bending mode. FAST, however, cannot yet model the shaft bending and, therefore, our results do not show any drivetrain whirl.  
+
+![](0ac7819b3ae0e7c2aacf262aac113c533dfc7ca943893a0ad6749a37484a84cc.jpg)  
+Figure 1.  Variation of in-vacuum modal frequencies with rotor speed (Coleman diagram).  
+
+As with the lag modes, the rotor exhibits flap modes that represent the cumulative effect of all the blades. Each mode, however, has a different physical significance and interacts differently with the fixed (tower-nacelle) system. The collective flap mode is a rotor mode in which all the blades flap, out of the rotor plane, in unison and with the same magnitude. One can visualize that in this mode, the rotor cones back and forth and interacts with the tower fore-aft mode. In the progressive flap mode, the rotor disk wobbles in the same direction as the rotor spin direction and at an angular speed higher than the rotor speed. In the regressive flap mode, the rotor disk wobbles in a direction opposite to that of the rotor spin direction and at an angular speed lower than the rotor speed.  
+
+In Fig. 1, we note that the nacelle yaw frequency, around $6.1\ \mathrm{Hz}$ , is hardly influenced by the rotor speed as expected. The variation of other frequencies is difficult to discern in this figure. Therefore, the results are replotted in Fig. 2. Note that the tower’s second side-to-side (Tower_SS2) modal frequency increases somewhat with rotor speed; this is due to its coupling with the second flap collective mode, which stiffens centrifugally as the rotor speed increases. Because of centrifugal stiffening, the frequencies of other collective modes (Flap1_coll and Lag1_coll) also increase with rotor speed. The frequencies of all progressive modes also increase with rotor speed. Note, however, that the progressive mode frequencies increase at rates much higher than those of the corresponding collective modes. All regressive mode frequencies decrease with rotor speed.  
+
+Readers familiar with the isolated rotor modes, wherein the tower-nacelle subsystem is fully rigid, will notice interesting differences. For an isolated rotor, progressive and regressive modal frequency curves are symmetric about the corresponding collective modal frequency curve; at any rotor speed, $\it{\Omega}_{\it{\Omega}}\it{\Omega}_{\it{\Omega}}$ the progressive mode frequency is $\omega\_c o l l\ +\mathcal{O},$ and the regressive mode frequency is ω_coll $-{\mathcal{Q}},$ where ω_coll is the collective mode frequency. However, this symmetry is lost for the full system because different rotor modes interact differently with the towernacelle subsystem. Even more interestingly, the Lag1_coll frequency curve is not even located between the Lag1_prog and Lag1_reg curves, normally observed for an isolated rotor. This is because the first rotor lag collective mode strongly interacts with the torsionally compliant drivetrain mode. The drivetrain mode frequency itself, however, is unaffected by the rotor speed. A more in-depth discussion of modal interaction amongst several subsystems is not in the scope of this paper.  
+
+![](f9a37716afcf03679f24d4f0c99575d23bc415c73f027d342f30895c406e2251.jpg)  
+Figure 2.  Variation of in-vacuum modal frequencies with rotor speed (land-based configuration).  
+
+![](b999f05672f6a4e076dbec24a08afbc2dbd6cc262957e5c25839277611cbef53.jpg)  
+Figure 3.  Variations of controlled blade-pitch angle and rotor speed with wind speed.  
+
+Now we consider normal operation in the presence of wind. The wind speed varies from $3~\mathrm{m/s}$ (cut-in speed) to $25~\mathrm{m/s}$ (cut-out speed). For a given wind speed, each blade’s pitch and rotor speed are set to the values depicted in Fig. 3. These values ensure optimum performance in region 2 and power regulation in region 3 [15]. Figs. 4 and 5 show the stability results.  
+
+![](b50bf56331e59666e941b360bf8e10e9a49c072c65e6e79e6f6deb2a6ff72fac.jpg)  
+Figure 4.  Variations of damped modal frequencies with wind speed (land-based configuration).  
+
+Figure 4 shows the variation of damped modal frequencies with wind speed. We make a few observations. Frequencies of the nacelle yaw, Flap2_coll, and Flap1_coll, and Flap1_prog modes decrease with wind speed, indicating a strong influence of aeroelastic couplings. A sudden change in the frequency trends near $11~\mathrm{m/s}$ wind speed is due to the fact that the rotor speed is held constant beyond this wind speed (see Fig. 3). Note that all tower modes are relatively insensitive to the wind speed.  
+
+![](ee9cd92dcccc22fd430a3bb04e98d449eff0153163a0ec5b4df597b11363aef9.jpg)  
+Figure 5.  Variations of modal damping ratios with wind speed (land-based configuration).  
+
+Figure 5 shows the variation of modal damping ratios with wind speed. All modes have positive damping, implying that they are all stable. As expected, all flap modes are highly damped. The damping of the rotor first flap modes (not shown in the plots) exceeding $50\%$ is somewhat high; in our experience with typical rotors, this damping is in the range $30\%{-}50\%$ . This suggests that the blade is quite light in comparison to its aerodynamic capability. Typical blades have a Lock number in the range 8–10 (a Lock number is a nondimensional parameter that expresses the aerodynamic lift capability of a blade in comparison to its weight [4]). The blade in this study has an effective Lock number of about 12. Damping of both the progressive and regressive lag modes increases rapidly as wind speed is increased from $3~\mathrm{m/s}$ to about $11~\mathrm{m/s}$ ; thereafter it decreases and becomes almost constant beyond the wind speed of $22\;\mathrm{m/s}$ . The drivetrain damping increases monotonically except near a wind speed of $12\;\mathrm{m}/\mathrm{s}$ where it shows a slight dip. The tower first side-to-side mode is the least damped showing maximum damping of about $0.65\%$ .  
+
+# VI.  Conclusions  
+
+A new multiblade coordinate transformation (MBC) scheme is presented that is applicable to both constant- and variable-speed turbines; it also allows dissimilar blades. The scheme covers control, disturbance, and output-related matrices, which are often required by control designers. Equations governing MBC are developed for both the firstand the second-order systems. Next, we showed that the common approach of first time-averaging the system matrices over the rotor azimuth and then applying eigenanalysis can lead to erroneous modal and stability predictions. Direct averaging eliminates all periodic terms that contribute to system dynamic; therefore causing erroneous results. While averaging itself is not a bad approach, it must follow MBC. Finally, we showed application of MBC to the modal and stability analysis of a 5-MW turbine.  
+
+Though the MBC formulation developed in this paper assumes a three-bladed turbine, the formulation is general enough to be readily extendable to a turbine with more than three blades. We plan to do so in the future.  
+
+# Acknowledgments  
+
+Thanks to Sandy Butterfield from NREL for motivating this effort, to Jason Jonkman for providing the 5-MW turbine data, and to Kathy O’Dell for editing this paper. This work was performed in support of the U.S. Department of Energy under contract number DE-AC36-83CH10093.  
+
+# References  
+
+1Miller, R. H., “Helicopter control and stability in hovering flight,” JAS , 15:8, August 1948.  
+
+2Coleman, R. P. and Feingold, A. M., “Theory of self-excited mechanical oscillations of helicopter rotors with hinged blades,” NACA Technical Report TR 1351, 1958. 3Hohenemser, K. H. and Yin, S. K., “Some applications of the method of multiblade coordinates,” Journal of the American Helicopter Society , Vol. 17, No. 3, July 1972. 4Johnson, W., Helicopter Theory. Princeton University Press, New Jersey, 1980. 5Bir, G., Chopra, I., et al. “University of Maryland Advanced Rotor Code (UMARC) Theory Manual,” Technical Report UM-AERO 94-18, Center for Rotorcraft Education and Research, University of Maryland, College Park, July 1994. 6Bir, G.S., Wright, A.D. and Butterfiled, S., “Stability Analysis of Variable-Speed Wind Turbines”, Proceedings of the 1997 ASME Wind Energy Symposium, Reno, January 6-9, 1997. 7Malcolm, D. J., “Modal Response of 3-Bladed Wind turbines,” Journal of Solar Engineering, Vol. 124, No. 4, Nov. 2002, pp. 372-377. 8Elliott, A.S.; McConville, J.B.  (1989). "Application of a General-Purpose Mechanical Systems Analysis Code to Rotorcraft Dynamics Problems."  Presented at the American Helicopter Society National Specialists' Meeting on Rotorcraft Dynamics, 1989. 9McCoy, T. J., “Wind turbine ADAMS model linearization including rotational and aerodynamic effects,” 2004 ASME Wind Energy Symposium, Jan. 2004, pp. 224-233. 10Hansen, M. H., “Improved modal dynamics of wind turbines to avoid stall-induced vibrations,” Wind Energy, Vol. 6, Issue 2, 2003, pp.179-195. 11Hansen, M. H., “Stability analysis of three-bladed turbines using an eigenvalue approach,” 2004 ASME Wind Energy Symposium, Jan. 2004, pp. 192-202. 12Riziotis, V. A., Voutsinas, S. G., Politis, E. S., and Chaviaropoulos, P. K., “Aeroelastic stability of wind turbines: The problem, the methods, and the issues,” Wind Energy, Vol. 7, Issue 4, 2004, pp.373-392. 13Jonkman, J. M. and Buhl, M. L., Jr., “FAST User’s Guide,” NREL/EL-500-29798, Golden, CO: National Renewable Energy Laboratory, October 2004. 14Bir, G. and Jonkman, J., “Aeroelastic instabilities of large offshore and onshore wind turbines,” Journal of Physics: Conference Series, 2007. 15 Jonkman, J.; Butterfield, S.; Musial, W.; and Scott, G., “Definition of a 5-MW Reference Wind Turbine for Offshore System Development,” NREL/TP-500-38060, Golden, CO: National Renewable Energy Laboratory, February 2007 (to be published).  
+
+16Bir, G. and Stol K., “Modal analysis of a teetered-rotor turbine using the Floquet approach,” Proceedings of the AIAA/ASME Wind Energy Symposium, Jan. 2000. 17 IEC 61400–1 Ed. 3, Wind Turbines – Part 1: Design Requirements, International Electrotechnical Commission (IEC), 2005. 18 IEC 61400–3, Wind Turbines – Part 3: Design Requirements for Offshore Wind Turbines, International Electrotechnical Commission (IEC), 2006 (to be published).  

线性化求解器/参考文献/Hansen-Improved modal dynamics of wind turbine/auto/Hansen-Improved modal dynamics of wind turbine.md

@@ -0,0 +1,470 @@
+# Improved Modal Dynamics of Wind Turbines to Avoid Stall-induced Vibrations  
+
+M. H. Hansen∗ Wind Energy Department, Risø National Laboratory, PO Box 49, DK-4000 Roskilde, Denmark  
+
+Key words: horizontal axis turbines; mechanical vibrations; modal analysis; optimization  
+
+Stall-induced edgewise blade vibrations have occasionally been observed on three-bladed wind turbines over the last decade. Experiments and numerical simulations have shown that these blade vibrations are related to certain vibration modes of the turbines. A recent experiment with a 600 kW turbine has shown that a backward whirling mode associated with edgewise blade vibrations is less aerodynamically damped than the corresponding forward whirling mode.Inthisarticlethemodeshapesoftheparticularturbineareanalysed, based on a simplified turbine model described in a multi-blade formulation. It is shown that the vibrations of the blades for the backward and forward edgewise whirling modes are different, which can explain the measured difference in aerodynamic damping. The modal dynamics of the entire turbine is important for stability assessments; blade-only analysis can be misleading. In some cases the modal dynamics may even be improved to avoid stall-induced vibrations. 
+
+过去十年中，在三叶风轮上偶有观察到失速诱发的摆振叶片振动。实验和数值模拟表明，这些叶片振动与机组的某些振动模态相关。最近对一台 600 千瓦风轮的实验表明，与摆振叶片振动相关的向后旋转模态的空气动力学阻尼小于相应的向前旋转模态。本文分析了特定风轮的模态形状，基于多叶片公式描述的简化风轮模型。结果表明，向后和向前摆振叶片振动的振动情况不同，这可以解释空气动力学阻尼的测量差异。整个风轮的模态动力学对于稳定性评估至关重要；仅对叶片进行分析可能会产生误导。在某些情况下，模态动力学甚至可以得到改善，以避免失速诱发的振动。
+
+Copyright $\circledcirc$ 2003 John Wiley & Sons, Ltd.  
+
+# Introduction  
+
+A theoretical modal analysis of a three-bladed wind turbine is presented, with the aim of explaining observed differences in aerodynamic damping of two apparently similar modes.1 Similar results have previously been published,2 but without the present details on the method and analysis.  
+一篇关于三叶风轮的理论模态分析被呈现，旨在解释两种看似相似的模态在气动阻尼上的观测差异。类似的结果已发表过，但缺乏本研究方法和分析的详细信息。
+
+Stall-induced vibrations of wind turbine blades have become an issue in the optimized design of stallregulated wind turbines with increasing rotor diameters. Until the 1990s, no problems with such vibrations had been observed, although work performed in the 1980s indicated that these problems could arise.3When one of the first examples of stall-induced vibrations was reported on a $500~\mathrm{kW}$ wind turbine,4 it led to further investigations which have provided a detailed understanding of the problem.1,5–8  
+
+From the studies of stall-induced vibrations, different types of solutions have been developed. The common solution for a wind turbine in operation is to change the stall characteristics of the blade by placing stall strips on the leading edge of the blade.8,9 When designing a new blade, this aerodynamic issue must of course be considered in the choice of aerofoil profiles. However, the structural dynamics of the blade is also important for the risk of stall-induced vibrations. The aerodynamic damping of a blade section operating in stall is highly dependent on the direction of blade vibration (relative to the rotor plane) associated with the primary bending modes.5,7,8 For most aerofoil profiles at high angles of attack it can be shown, using quasi-steady aerodynamics (see Appendix), that vibrations in the rotor plane (edgewise) are less damped by the air flow than those out of the rotor plane (flapwise). Another important role of the structural dynamics is the structural damping of the blade. For example, if the first edgewise mode of a blade has a negative aerodynamic damping at some wind speeds, then it may be possible to avoid an instability by adding sufficient structural damping to the blade, or a distinct vibration damper.10 Thus, besides the aerodynamic issue, it is necessary to consider the mode shapes and their structural damping when designing a new blade.  
+
+Another issue of stall-induced vibrations is the interaction of the blades with the remaining wind turbine. A numerical simulation of a turbine using the aeroelastic code HAWC11 showed that, by stiffening the support of the rotor (shaft and nacelle), undamped edgewise blade vibrations could be avoided.7 This observation may be related to an observation made by Thomsen et al.1 They have developed an experimental method for estimating the total damping of a wind turbine during operation. Figure 1 shows the results of their experiments with a stall-controlled Bonus $600~\mathrm{kW}$ wind turbine. The observation is that the forward rotor whirling mode associated with edgewise blade vibrations is more damped than the backward rotor whirling mode. The structural damping of these two modes is assumed to be the same, because their natural frequencies and mode shapes are almost identical; the difference in total damping is therefore a difference in aerodynamic damping.  
+
+In this article it is suggested that the dependence of aerodynamic damping on whirl direction and rotor support stiffness can be explained by analysing the turbine mode shapes. It is shown how a theoretical modal analysis can be performed for a rotating turbine based on the so-called multi-blade co-ordinate transformation.12,13 For the particular $600~\mathrm{kW}$ turbine it is shown that the blades vibrate more out of the rotor plane in the forward edgewise whirling mode than in the backward edgewise whirling mode. This observation can explain the measured difference in aerodynamic damping for the two whirling modes, because the direction of blade vibration with respect to the rotor plane is important for the aerodynamic damping in stall. It is suggested that, using the modal analysis, the design of wind turbines can be tailored to increase the blade vibration out of the rotor plane for all modes, and thereby increase the aerodynamic damping of stall-regulated turbines.  
+
+The next section describes the model and method for modal analysis of rotating turbines. The modal analysis of the $600~\mathrm{kW}$ turbine is then presented, followed by a section on the measured differences in aerodynamic damping and the optimization of its dynamical behaviour.  
+
+随着风轮直径不断增大，叶片失速诱导振动已成为优化失速调控风力发电机设计的难题。直到 1990 年代之前，人们并未观察到此类振动问题，尽管 1980 年代的研究表明这些问题可能会出现。<sup>3</sup> 当第一个失速诱导振动案例在 500 千瓦的风力发电机上被报告时，<sup>4</sup> 它引发了进一步的调查，这些调查提供了对该问题的详细理解。<sup>1,5–8</sup>
+
+针对失速诱导振动，已经开发出多种解决方案。风力发电机运行时的常见解决方案是在叶片前缘放置失速条，从而改变叶片的失速特性。<sup>8,9</sup> 在设计新叶片时，必须在气翼剖面选择中考虑这一空气动力学问题。然而，叶片的结构动力学也对失速诱导振动的风险至关重要。叶片在失速状态下运行的叶片段的空气动力学阻尼高度依赖于叶片振动方向（相对于风轮平面）与主要弯曲模态之间的关系。<sup>5,7,8</sup> 对于大多数气翼剖面在较高攻角时，使用准稳态空气动力学（见附录）可以证明，风轮平面内的振动（摆振）比风轮平面外的振动（挥舞）受气流阻尼较小。结构动力学的另一个重要作用是叶片的结构阻尼。例如，如果叶片的第一个摆振模态在某些风速下具有负的空气动力学阻尼，那么可以通过向叶片添加足够的结构阻尼或一个独立的振动阻尼器来避免不稳定性。<sup>10</sup> 因此，除了空气动力学问题，在设计新叶片时，必须考虑模态形状及其结构阻尼。
+
+失速诱导振动的另一个问题是叶片与剩余风力发电机之间的相互作用。使用气弹耦合代码 HAWC11 的涡轮数值模拟表明，通过加固风轮（轴和舱室）的支撑，可以避免未阻尼的摆振叶片振动。<sup>7</sup> 这一观察结果可能与 Thomsen 等人所做的观察相关。<sup>1</sup> 他们开发了一种实验方法来估计风力发电机在运行期间的总阻尼。图 1 显示了他们对失速调控 Bonus 600 千瓦风力发电机进行的实验结果。观察结果是，与摆振叶片振动相关的正向风轮回转模态比反向风轮回转模态更受阻尼。假设这两个模态的结构阻尼相同，因为它们的固有频率和模态形状几乎相同；因此，总阻尼的差异是空气动力学阻尼的差异。
+
+本文建议，空气动力学阻尼对回转方向和风轮支撑刚度的依赖性可以通过分析涡轮模态形状来解释。展示了如何基于所谓的“多叶片坐标变换”<sup>12,13</sup> 对旋转涡轮进行理论模态分析。对于特定的 600 千瓦涡轮，表明在正向摆振回转模态中，叶片比在反向摆振回转模态中更倾向于风轮平面外振动。这一观察结果可以解释两种回转模态的空气动力学阻尼差异的测量结果，因为叶片相对于风轮平面的振动方向对失速状态下的空气动力学阻尼至关重要。建议使用模态分析，可以调整风力发电机设计，以增加所有模态的风轮平面外叶片振动，从而增加失速调控涡轮的空气动力学阻尼。
+
+下一节描述了旋转风机模态分析的模型和方法。然后，对 600 千瓦涡轮进行模态分析，随后是关于测量到的空气动力学阻尼差异和其动力学行为优化的部分。
+
+
+# Model and Method  
+
+A linear model of a three-bladed wind turbine is needed for the modal analysis. It must qualitatively describe the dynamics of the turbine up to the second tilt/yaw modes in order to study the phenomena mentioned above. Blade torsion is not included in the analysis, because the corresponding natural frequencies are several times higher than the dynamics of interest.  
+为了模态分析，需要建立一个三叶风电机组的线性模型。为了研究上述现象，它必须定性描述风电机组直至第二倾覆/偏航模态的动力学。由于其固有频率远高于感兴趣的动力学，因此未在分析中包含叶片扭转。
+
+
+![](146a343862b1487203627349c8c528f86914d6bbb025069461ed2230d451d999.jpg)  
+Figure 1. Total damping of the forward and backward edgewise whirling modes of the Bonus 600 kW estimated by Thomsen et al.1  
+
+![](c54c4be6a0b3a29c3948930025463b595b9a6124393996423d58702f62526d82.jpg)  
+Figure 2. Simplified structural model of a wind turbine  
+
+# Equations of Motion  
+
+A schema of the structural model is shown in Figure 2. The motion of the nacelle and tower is described by seven degrees of freedom (DOFs). The nacelle can translate in the two horizontal directions, described in the ground fixed frame $(X,Y,Z)$ by $u_{x}^{\mathrm{t}}$ (lateral) and $u_{y}^{\mathrm{t}}$ (longitudinal), and can tilt and yaw about the tower top, described by the angles $\theta_{x}^{\mathrm{t}}$ and $\theta_{z}^{\mathrm{t}}$ respectively. The rotor can rotate out of its plane owing to bending of the shaft, described by the angles $\theta_{x}^{\mathrm{s}}$ and $\theta_{z}^{\mathrm{s}}$ defined in the corotating frame of blade number 1. Finally, the torsion of the shaft and drive-train is described by the angle $\theta_{y}^{\mathrm{s}}$ .  
+
+The motions of the blades are described by a modal expansion defined in their own co-rotating frames $(x,\,y,\,z)$ , where the $z_{i}$ -axis is the blade axis and the $y$ -axis at rest coincides with the $Y$ -axis. The in- and out-of-rotor-plane motions of blade number $i$ are given by  
+
+图2所示为结构模型的示意图。机舱和塔架的运动由七个自由度（DOF）描述。机舱可以在两个水平方向上平移，分别在地面固定坐标系 $(X,Y,Z)$ 中描述为 $u_{x}^{\mathrm{t}}$（X）和 $u_{y}^{\mathrm{t}}$（Y），并且可以绕塔顶倾覆和偏航，分别由角度 $\theta_{x}^{\mathrm{t}}$ 和 $\theta_{z}^{\mathrm{t}}$ 描述。风轮由于轴的弯曲而可以偏离其平面，分别由角度 $\theta_{x}^{\mathrm{s}}$ 和 $\theta_{z}^{\mathrm{s}}$ 在叶片编号1的共转坐标系中定义。最后，轴和驱动系的扭转由角度 $\theta_{y}^{\mathrm{s}}$ 描述。
+
+叶片的运动由在其自身共转坐标系 $(x,\,y,\,z)$ 中定义的模态展开来描述，其中 $z_{i}$ 轴为叶片轴，静止时的 $y$ 轴与 $Y$ 轴重合。叶片编号 $i$ 的进出风轮平面运动由给出。
+
+$$
+u_{x}^{\mathrm{B}_{i}}(t,z)=\sum_{n=1}^{N}q_{i,n}(t)\Phi_{n}^{x}(z),\qquad u_{y}^{\mathrm{B}_{i}}(t,z)=\sum_{n=1}^{N}q_{i,n}(t)\Phi_{n}^{y}(z)
+$$  
+
+where $q_{i,n}$ are modal blade $c o$ -ordinates describing the contents of mode number $n$ in the motion of blade number $i$ , and $\Phi_{n}^{x}$ and $\Phi_{n}^{y}$ are the mode shape functions. These shape functions can be obtained numerically or experimentally.  
+
+For the derivation of the inertia forces it is assumed that the centre of mass (CM) for the undeformed blade lies along the $z_{i}$ -axis. A vector in the blade frame $(x,\,y,z)$ from the rotor centre to the CM at radius $z$ is therefore given by  
+其中，$q_{i,n}$ 为模态叶片坐标，描述了第 $i$ 叶片运动中第 $n$ 模态的内容，而 $\Phi_{n}^{x}$ 和 $\Phi_{n}^{y}$ 为模态形函数。这些形函数可以通过数值计算或实验获得。
+
+在推导惯性力时，假设未变形叶片的质量中心（CM）位于 $z_{i}$ 轴上。因此，从机组中心到半径为 $z$ 的质量中心的叶片坐标系 $(x,\,y,z)$ 中的向量表示为：
+$$
+\mathbf{r}_{i}^{\mathrm{B}}=\sum_{n=1}^{N}\{\Phi_{n}^{x}(z),\,\Phi_{n}^{y}(z),z\}^{\mathrm{T}}q_{i,n}(t)
+$$  
+
+In the ground fixed frame $(X,Y,Z)$ the motion of the CM at radius $z$ on blade number $i$ becomes  
+在固定地面坐标系 $(X,Y,Z)$ 中，半径为 $z$，叶片编号为 $i$ 的质心运动变为：
+
+$$
+\mathbf{r}_{i}=\mathbf{r}_{\mathrm{{t}}}+\mathbf{T}_{\mathrm{{t}}}[\mathbf{r}_{\mathrm{{ts}}}+\mathbf{T}_{\mathrm{{a}}}\mathbf{T}_{\mathrm{{s}}}(\mathbf{r}_{\mathrm{{sh}}}+\mathbf{T}_{\mathrm{{B}}_{i}}\mathbf{r}_{i}^{\mathrm{{B}}})]
+$$  
+
+where ${\bf r}_{\mathrm{t}}=\{u_{x}^{\mathrm{t}},u_{y}^{\mathrm{t}},0\}^{\mathrm{T}}$ is a vector from the origin to the tower top, ${\bf r}_{\mathrm{ts}}=\{0,-L_{\mathrm{n}},0\}^{\mathrm{T}}$ is a vector from the tower top to the bending point on the shaft in nacelle co-ordinates, and $\mathbf{r}_{\mathrm{sh}}=\{0,-L_{\mathrm{s}},0\}^{\mathrm{T}}$ is a vector from this point to the rotor centre in rotor co-ordinates. The transformation matrices ${\mathbf{T}}_{\mathrm{a}}$ and ${\bf{T}}_{{\bf{B}}_{i}}$ handling the azimuth rotations are written as  
+
+其中 ${\bf r}_{\mathrm{t}}=\{u_{x}^{\mathrm{t}},u_{y}^{\mathrm{t}},0\}^{\mathrm{T}}$ 是从原点到塔顶的矢量，${\bf r}_{\mathrm{ts}}=\{0,-L_{\mathrm{n}},0\}^{\mathrm{T}}$ 是在机舱坐标系中，从塔顶到曲柄轴弯曲点的矢量，而 $\mathbf{r}_{\mathrm{sh}}=\{0,-L_{\mathrm{s}},0\}^{\mathrm{T}}$ 是从该点到风轮中心的矢量，坐标系为风轮坐标系。处理方位角旋转的变换矩阵 ${\mathbf{T}}_{\mathrm{a}}$ 和 ${\bf{T}}_{{\bf{B}}_{i}}$ 如下：
+$$
+\begin{array}{r l}&{\mathbf{T}_{\mathbf{B}_{i}}=\left[\begin{array}{c c c}{\cos\left[\frac{2}{3}\pi(i-1)\right]}&{0}&{\sin\left[\frac{2}{3}\pi(i-1)\right]}\\ {0}&{1}&{0}\\ {-\sin\left[\frac{2}{3}\pi(i-1)\right]}&{0}&{\cos\left[\frac{2}{3}\pi(i-1)\right]}\end{array}\right]}\\ &{\mathbf{T}_{\mathbf{a}}=\left[\begin{array}{c c c}{\cos(\Omega t+\theta_{\mathbf{y}}^{\mathrm{s}})}&{0}&{\sin(\Omega t+\theta_{\mathbf{y}}^{\mathrm{s}})}\\ {0}&{1}&{0}\\ {-\sin(\Omega t+\theta_{\mathbf{y}}^{\mathrm{s}})}&{0}&{\cos(\Omega t+\theta_{\mathbf{y}}^{\mathrm{s}})}\end{array}\right]}\end{array}
+$$  
+
+where $\Omega$ is the rotation speed. It is assumed that all rotations of the tower top and hub are small $(\theta_{x}^{\mathrm{t}},\theta_{z}^{\mathrm{t}},\theta_{x}^{\mathrm{s}},\theta_{z}^{\mathrm{s}}<<1)$ , whereby the remaining transformation matrices of equation (3) become  
+其中 $\Omega$ 为转速。假设塔顶和轮毂的所有旋转均较小 $(\theta_{x}^{\mathrm{t}},\theta_{z}^{\mathrm{t}},\theta_{x}^{\mathrm{s}},\theta_{z}^{\mathrm{s}}<<1)$ ，因此方程(3)中剩余的变换矩阵变为
+$$
+\begin{array}{r}{\mathbf{T}_{\mathrm{s}}=\left[\begin{array}{l l l}{1}&{-\theta_{z}^{\mathrm{s}}}&{0}\\ {\theta_{z}^{\mathrm{s}}}&{1}&{-\theta_{x}^{\mathrm{s}}}\\ {0}&{\theta_{x}^{\mathrm{s}}}&{1}\end{array}\right],\ \mathbf{T}_{\mathrm{t}}=\left[\begin{array}{l l l}{1}&{-\theta_{z}^{\mathrm{t}}}&{0}\\ {\theta_{z}^{\mathrm{t}}}&{1}&{-\theta_{x}^{\mathrm{t}}}\\ {0}&{\theta_{x}^{\mathrm{t}}}&{1}\end{array}\right]}\end{array}
+$$  
+
+where ${\bf{T}}_{\mathrm{{s}}}$ handles the transformation from rotor to nacelle co-ordinates due to shaft bending, and ${\bf T}_{\mathrm{t}}$ handles the transformation from nacelle co-ordinates to the ground fixed frame due to tilt and yaw.  
+
+The total kinetic energy of the wind turbine in this model can now be derived as  
+其中 ${\bf{T}}_{\mathrm{{s}}}$ 处理由于主轴弯曲引起的风轮到机舱坐标系的变换，而 ${\bf T}_{\mathrm{t}}$ 处理由于倾覆和偏航引起的机舱坐标系到地面固定坐标系的变换。
+
+在此模型中，风电机组的总体动能现在可以推导如下：
+$$
+T=\frac{1}{2}\left(M(\dot{u}_{x}^{t^{2}}+\dot{u}_{y}^{t^{2}})+I_{x}\dot{\theta}_{x}^{t^{2}}+I_{z}\dot{\theta}_{z}^{t^{2}}+I_{y}\dot{\theta}_{y}^{s^{2}}+\sum_{i=1}^{3}\int_{0}^{R}m(z)|\dot{\bf r}_{i}|^{2}{\bf d}z\right)
+$$  
+
+where $\dot{\left(\,\right)}\equiv\partial/\partial t$ denotes the time derivative, $M$ is the total mass of the nacelle (and a part of the tower), $I_{x}$ and $I_{z}$ are the rotational inertias of the nacelle about the tilt and yaw axes respectively, $I_{y}$ is the rotational inertia of the shaft and drive-train, $R$ is the length of the blades including the hub and root extension, and $m(z)$ is the mass per unit length of the entire blade–root–hub assembly.  
+
+It is assumed that there is an elastic stiffness coupling between the nacelle tilt $\theta_{x}^{\mathrm{t}}$ and the longitudinal tower motion $u_{y}^{\mathrm{t}}$ , whereby the potential energy is defined as  
+
+其中 $\dot{\left(\,\right)}\equiv\partial/\partial t$ 表示时间导数，$M$ 为机舱（及塔架的一部分）的总质量，$I_{x}$ 和 $I_{z}$ 分别为机舱绕tilt轴和偏航轴的转动惯量，$I_{y}$ 为轴和传动系的转动惯量，$R$ 为叶片（包括轮毂和根部延伸）的长度，$m(z)$ 为整个叶片-根部-轮毂组件的单位长度质量。
+
+假设机舱tilt角 $\theta_{x}^{\mathrm{t}}$ 和纵向塔架运动 $u_{y}^{\mathrm{t}}$ 之间存在弹性刚度耦合，势能定义如下：
+$$
+\begin{array}{r l r}{{V=\frac{1}{2}\left[k_{x}u_{x}^{\mathrm{t}^{2}}+k_{y}u_{y}^{\mathrm{t}^{2}}+g_{x y}u_{y}^{\mathrm{t}}\theta_{x}^{\mathrm{t}}+G_{x}\theta_{x}^{\mathrm{t}^{2}}+G_{z}\theta_{z}^{\mathrm{t}^{2}}\right.}}\\{+K_{s}(\theta_{x}^{s^{2}}+\theta_{z}^{s^{2}})+G_{y}\theta_{y}^{s^{2}}+\sum_{i=1}^{3}\left(\sum_{n=1}^{N}k_{n}q_{i,n}^{2}+V_{\mathrm{c},i}\right)]}\end{array}
+$$  
+
+where $k_{x}$ and $k_{y}$ are the translational stiffnesses of the tower at its top, $g_{x y}$ is the coupling stiffness between nacelle tilt and longitudinal tower motion, $G_{x}$ and $G_{z}$ are the rotational stiffnesses of the tower at its top, $K_{\mathrm{s}}$ is the bending stiffness of the shaft, and $G_{y}$ is the torsional stiffness of the shaft and drive-train. The modal stiffness $k_{n}$ is defined as $k_{n}\equiv\omega_{n}^{2}M_{n}$ , where $\omega_{n}$ is the natural frequency of blade mode number $n$ , and $M_{n}$ is the modal mass, $\begin{array}{r}{M_{n}\equiv\int_{0}^{R}m(z)(\Phi_{n}^{x^{2}}+\Phi_{n}^{y^{2}})\mathbf{d}z}\end{array}$ .  
+
+The potential energy term $V_{\mathrm{c},i}$ is added for modelling of the centrifugal stiffening effect. The tensile force in the blade due to centrifugal body forces is $\begin{array}{r}{\Omega^{2}\int_{z}^{R}m(\zeta)\zeta\mathrm{d}\zeta}\end{array}$ , and the additional potential energy due to the centrifugal body forces can be approximated by  
+
+其中，$k_{x}$ 和 $k_{y}$ 分别为塔顶的平动刚度，$g_{x y}$ 为机舱倾斜与纵向塔架运动之间的耦合刚度，$G_{x}$ 和 $G_{z}$ 分别为塔顶的转动刚度，$K_{\mathrm{s}}$ 为轴的弯曲刚度，$G_{y}$ 为轴和驱动系统的扭转刚度。模态刚度 $k_{n}$ 定义为 $k_{n}\equiv\omega_{n}^{2}M_{n}$ ，其中 $\omega_{n}$ 为叶片模态数 $n$ 的固有频率， $M_{n}$ 为模态质量，$\begin{array}{r}{M_{n}\equiv\int_{0}^{R}m(z)(\Phi_{n}^{x^{2}}+\Phi_{n}^{y^{2}})\mathbf{d}z}\end{array}$ 。
+
+为了模拟离心加固效应，添加了势能项 $V_{\mathrm{c},i}$。由于离心惯性力，叶片中的拉力为 $\begin{array}{r}{\Omega^{2}\int_{z}^{R}m(\zeta)\zeta\mathrm{d}\zeta}\end{array}$ ，由于离心惯性力产生的附加势能可以近似为
+$$
+V_{\mathrm{c},i}=\frac{1}{2}\Omega^{2}\int_{0}^{R}\left[\left(\sum_{n=1}^{N}q_{i,n}\frac{\mathrm{d}\Phi_{n}^{x}}{\mathrm{d}z}\right)^{2}+\left(\sum_{n=1}^{N}q_{i,n}\frac{\mathrm{d}\Phi_{n}^{y}}{\mathrm{d}z}\right)^{2}\right]\int_{z}^{R}m(\zeta)\zeta\mathrm{d}\zeta\mathrm{d}z
+$$  
+where it is assumed that rotations of the blade cross-sections due to bending are small. 
+The blade and tower/nacelle co-ordinates are collected in a state vector  
+其中假设因弯曲引起的叶片横截面旋转较小。
+叶片和塔架/机舱坐标被收集到一个状态向量中。
+
+$$
+\mathbf{x}=\{q_{1,1},\dotsc,q_{1,N},q_{2,1},\dotsc,q_{2,N},q_{3,1},\dotsc,q_{3,N},\theta_{x}^{\mathrm{s}},\theta_{z}^{\mathrm{s}},u_{x}^{\mathrm{t}},u_{y}^{\mathrm{t}},\theta_{x}^{\mathrm{t}},\theta_{z}^{\mathrm{s}},\theta_{y}^{\mathrm{s}}\}^{\mathrm{T}}
+$$  
+
+and Lagrange’s equations are derived from the Lagrangian $L=T-V$ . After linearization about the undeformed state by assuming small amplitudes of the co-ordinates, $x_{j}<<1$ , the equations of motion can be written in the form  
+并且拉格朗日方程可由拉格朗日量 $L=T-V$ 推导得出。在假设坐标变量 $x_{j}$ 具有小幅值，即 $x_{j}<<1$，并对未变形状态进行线性化后，运动方程可以写成如下形式：
+$$
+\mathbf{M}(t){\ddot{\mathbf{x}}}+\mathbf{C}(t){\dot{\mathbf{x}}}+\mathbf{K}(t)\mathbf{x}=\mathbf{0}
+$$  
+
+Where M, C and $\mathbf{K}$ are the mass, gyroscopic and stiffness matrices respectively. All three matrices are time-dependent, containing harmonic terms with period $T=2\pi/\Omega$ .  
+
+Note that structural damping is neglected in this simplified wind turbine model. The lower-order modes of interest to this study are low damped, thus structural damping has practically no effect on their mode shapes. For studies of resonance phenomena and stability it would be necessary to include the structural damping.14  
+
+The equations of motion (10) have a form that is general for linear structural models of wind turbines. The periodic terms arising from the azimuthal rotation prevent a direct derivation of the natural frequencies and mode shapes from an eigenvalue problem. Still, the mathematical approach of using Floquet theory should be avoided, because there is a more physically comprehensive way to set up an eigenvalue problem. 
+
+其中 M、C 和 $\mathbf{K}$ 分别为质量、陀螺惯性力和刚度矩阵。这三个矩阵均随时间变化，包含周期为 $T=2\pi/\Omega$ 的谐波项。
+在简化的风电机组模型中，忽略了结构阻尼。本研究感兴趣的低阶模态阻尼较小，因此结构阻尼实际上对它们的模态形状没有影响。若要研究共振现象和稳定性，则需要包含结构阻尼。<sup>14</sup>
+运动方程 (10) 具有适用于风电机组线性结构模型的通用形式。由于方位角旋转产生的周期项，无法从特征值问题直接推导出固有频率和模态。然而，仍应避免使用Floquet理论进行数学处理，因为存在一种更具物理意义的方式来设置特征值问题。 
+
+# Multi-blade Co-ordinate Transformation  
+
+The use of multi-blade co-ordinates, or Fourier co-ordinates, for bladed rotors12,13 is a method to describe the motions of individual blades in the same co-ordinate system as the structure supporting the rotor, whereby the periodic terms in the governing equations are eliminated. The fundamental assumption of this method is that the rotor is isotropic, i.e. all blades are identical, identically pitched and symmetrically mounted on the hub.  
+
+In the wind turbine model given by equation (10), the motion of the tower/nacelle is described in the ground fixed frame, while the shaft bending is described in rotating shaft co-ordinates, and the motions of the three blades are described in their own co-rotating frames. The shaft bending can easily be described in the ground fixed frame by a co-ordinate transformation with respect to the azimuth angle. The transformation of the blade co-ordinates into the ground fixed frame can be done by the multi-blade co-ordinate transformation, which is defined by  
+
+使用多叶片坐标，或傅里叶坐标，来描述叶片风轮的运动<sup>12,13</sup>是一种方法，它将单个叶片的运动描述在支撑风轮的结构相同的坐标系下，从而消除了控制方程中的周期项。该方法的根本假设是风轮是各向同性的，即所有叶片都相同，变桨角度相同，并且对称地安装在轮毂上。
+
+在由公式(10)给出的风电机组模型中，塔架/机舱的运动在地面固定坐标系中描述，主轴弯曲在旋转主轴坐标系中描述，而三个叶片的运动则在各自的共旋转坐标系中描述。通过相对于方位角的坐标变换，主轴弯曲可以很容易地在地面固定坐标系中描述。叶片坐标变换到地面固定坐标系可以通过多叶片坐标变换来实现，该变换由以下定义：
+
+$$
+q_{i,n}=a_{0,n}+a_{1,n}\cos\psi_{i}+b_{1,n}\sin\psi_{i}
+$$  
+
+where $\begin{array}{r}{\psi_{i}=\Omega t+\frac{2}{3}\pi(i-1)}\end{array}$ is the azimuth angle to blade number $i$ (without shaft torsion) measured from the $Z$ -axis.  
+其中 $\begin{array}{r}{\psi_{i}=\Omega t+\frac{2}{3}\pi(i-1)}\end{array}$ 是以 $Z$ 轴为基准，测得的叶片编号 $i$ 的方位角（无主轴扭转）。
+
+The three multi-blade co-ordinates $a_{0,n},\,a_{1,n}$ and $b_{1,n}$ replace the three blade co-ordinates $q_{1,n},\;q_{2,n}$ and $q_{3,n}$ . To see that the multi-blade co-ordinates describe the rotor motion in the ground fixed frame, assume that blade mode $n$ is the first flapwise bending mode. Then $a_{0,n}$ describes a simultaneous flapwise deflection of all three blades, while $a_{1,n}$ and $b_{1,n}$ describe tilt and yaw motions respectively.  
+
+三片叶片坐标 $a_{0,n},\,a_{1,n}$ 和 $b_{1,n}$ 取代了三片叶片坐标 $q_{1,n},\;q_{2,n}$ 和 $q_{3,n}$ 。为了说明多片叶片坐标描述了风轮在地面固定参考系中的运动，假设叶片模态 $n$ 是第一挥舞摆振模态。那么 $a_{0,n}$ 描述了所有三片叶片的同步挥舞变形，而 $a_{1,n}$ 和 $b_{1,n}$ 分别描述了倾覆和偏航运动。
+
+The multi-blade transformation of the state vector $\mathbf{x}$ can be represented by  
+
+$$
+\mathbf{x}=\mathbf{B}(t)\mathbf{z}
+$$  
+
+where the state vector $\mathbf{z}$ containing the multi-blade co-ordinates is  
+
+$$
+{\mathbf z}=\{a_{0,1},\ldots,a_{0,N},a_{1,1},\ldots,a_{1,N},b_{1,1},\ldots,b_{1,N},\theta_{X}^{\mathrm{s}},\theta_{Z}^{\mathrm{s}},u_{x}^{\mathrm{t}},u_{y}^{\mathrm{t}},\theta_{x}^{\mathrm{t}},\theta_{z}^{\mathrm{s}},\theta_{y}^{\mathrm{s}}\}^{\mathrm{T}}
+$$  
+
+and the transformation matrix is  
+
+$$
+\mathbf{B}\equiv\left[\begin{array}{c c c c c c}{\mathbf{I}_{N}}&{\mathbf{I}_{N}\cos\psi_{1}}&{\mathbf{I}_{N}\sin\psi_{1}}&{\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}\\ {\mathbf{I}_{N}}&{\mathbf{I}_{N}\cos\psi_{2}}&{\mathbf{I}_{N}\sin\psi_{2}}&{\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}\\ {\mathbf{I}_{N}}&{\mathbf{I}_{N}\cos\psi_{3}}&{\mathbf{I}_{N}\sin\psi_{3}}&{\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}\\ {\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}&{\cos\Omega t}&{-\sin\Omega t}&{\mathbf{0}}\\ {\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}&{\sin\Omega t}&{\cos\Omega t}&{\mathbf{0}}\\ {\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}&{\mathbf{I}_{S}.}\end{array}\right]
+$$  
+
+where ${\mathbf{I}}_{N}$ is the identity matrix of size $N\times N$ . The matrix $\mathbf{B}$ has the following properties:  
+
+$$
+\mathbf{B}^{-1}{\dot{\mathbf{B}}}={\left[\begin{array}{l l l l l l}{\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}\\ {\mathbf{0}}&{\mathbf{0}}&{{\boldsymbol{\Omega}}\mathbf{I}_{N}}&{\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}\\ {\mathbf{0}}&{-{\boldsymbol{\Omega}}\mathbf{I}_{N}}&{\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}\\ {\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}&{0}&{-{\boldsymbol{\Omega}}}&{\mathbf{0}}\\ {\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}&{{\boldsymbol{\Omega}}}&{0}&{\mathbf{0}}\\ {\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}&{\mathbf{0}}\end{array}\right]}\equiv\mathbf{R}
+$$  
+
+$$
+{\bf B}^{-1}=\mu{\bf B}^{\mathrm{T}},\quad\pmb{\mu}=\left[\begin{array}{c c c c}{\frac{1}{3}{\bf I}_{N}}&{{\bf0}}&{{\bf0}}&{{\bf0}}\\ {{\bf0}}&{\frac{2}{3}{\bf I}_{N}}&{{\bf0}}&{{\bf0}}\\ {{\bf0}}&{{\bf0}}&{\frac{1}{3}{\bf I}_{N}}&{{\bf0}}\\ {{\bf0}}&{{\bf0}}&{{\bf0}}&{{\bf I}_{7}}\end{array}\right]
+$$  
+
+which shows that $\dot{\mathbf{B}}=\mathbf{B}\mathbf{R}$ and $\ddot{\mathbf{B}}=\mathbf{B}\mathbf{R}^{2}$ . A transformation of equation (10) into the multi-blade co-ordinates (12) yields  
+
+$$
+\mathbf{M}_{\mathrm{B}}\ddot{\mathbf{z}}+(2\mathbf{M}_{\mathrm{B}}\mathbf{R}+\mathbf{C}_{\mathrm{B}})\dot{\mathbf{z}}+(\mathbf{M}_{\mathrm{B}}\mathbf{R}^{2}+\mathbf{C}_{\mathrm{B}}\mathbf{R}+\mathbf{K}_{\mathrm{B}})\mathbf{z}=0
+$$  
+
+where $\mathbf{M}_{\mathrm{B}}\equiv\mu\mathbf{B}^{\mathrm{T}}\mathbf{M}\mathbf{B}$ , $\mathbf{C}_{\mathrm{B}}\equiv\mu\mathbf{B}^{\mathrm{T}}\mathbf{C}\mathbf{B}$ and $\mathbf{K}_{\mathrm{B}}\equiv\mu\mathbf{B}^{\mathrm{T}}\mathbf{K}\mathbf{B}$ are the transformed system matrices. If the rotor is isotropic, these matrices are constant, i.e. the transformed equations of motion in the multi-blade co-ordinates (16) contain no periodic terms. It is therefore possible to define an eigenvalue problem that gives the natural frequencies and mode shapes of the wind turbine with a rotating rotor.  
+其中 $\mathbf{M}_{\mathrm{B}}\equiv\mu\mathbf{B}^{\mathrm{T}}\mathbf{M}\mathbf{B}$ , $\mathbf{C}_{\mathrm{B}}\equiv\mu\mathbf{B}^{\mathrm{T}}\mathbf{C}\mathbf{B}$ 和 $\mathbf{K}_{\mathrm{B}}\equiv\mu\mathbf{B}^{\mathrm{T}}\mathbf{K}\mathbf{B}$ 是变换后的系统矩阵。如果风轮是各向同性的，这些矩阵是常数，即在多叶片坐标系（16）中变换后的运动方程不包含周期项。因此，可以定义一个特征值问题，以获得具有旋转风轮的风电机组的固有频率和模态。
+# Eigenvalue Problem and Solution  
+
+A substitution of the solution $\mathbf{z}=\mathbf{z}_{k}\mathbf{e}^{\lambda_{k}t}$ into (16) yields an eigenvalue problem that determines the complex eigenvalues $\lambda_{k}\equiv\sigma_{k}+\mathrm{i}\omega_{k}$ $(\mathbf{i}\equiv{\sqrt{-1}})$ and complex eigenvectors $\mathbf{Z}_{k}$ of mode number $k$ .  
+Inserting the solution $\mathbf{z}=\mathbf{z}_{k}\mathbf{e}^{\lambda_{k}t}$ into (12), the blade co-ordinate $q_{i,n}$ for the modal content of blade mode number $n$ in the motion of blade number $i$ becomes  
+
+将解 $\mathbf{z}=\mathbf{z}_{k}\mathbf{e}^{\lambda_{k}t}$ 代入 (16) 可得到一个特征值问题，该问题确定了模态数 $k$ 的复特征值 $\lambda_{k}\equiv\sigma_{k}+\mathrm{i}\omega_{k}$ $(\mathbf{i}\equiv{\sqrt{-1}})$ 和复特征向量 $\mathbf{Z}_{k}$。
+将解 $\mathbf{z}=\mathbf{z}_{k}\mathbf{e}^{\lambda_{k}t}$ 代入 (12)，风轮叶片（blade）的坐标 $q_{i,n}$，用于描述第 $i$ 个叶片的第 $n$ 个模态（mode）运动的模态内容，变为：
+
+$$
+\begin{array}{r l}&{q_{i,n,k}=A_{n,k}^{0}\mathrm{e}^{\sigma_{k}t}\cos(\omega_{k}t+\phi_{\beta,k}^{0})}\\ &{\quad\quad\quad+\frac{1}{2}A_{n,k}^{\mathrm{BW}}\mathrm{e}^{\sigma_{k}t}\cos[(\omega_{k}+\Omega)t+\frac{2}{3}\pi(i-1)+\phi_{n,k}^{\mathrm{BW}}]}\\ &{\quad\quad\quad+\frac{1}{2}A_{n,k}^{\mathrm{FW}}\mathrm{e}^{\sigma_{k}t}\cos[(\omega_{k}-\Omega)t-\frac{2}{3}\pi(i-1)+\phi_{n,k}^{\mathrm{FW}}]}\end{array}
+$$  
+
+where the amplitudes and phases are determined by the complex eigenvectors  
+其中，幅值和相位由复特征向量决定。
+$$
+\begin{array}{r l}&{A_{n,k}^{0}=\frac{1}{2}\sqrt{\mathrm{Re}\{a_{0,n,k}\}^{2}+\mathrm{Im}\{a_{0,n,k}\}^{2}}}\\ &{A_{n,k}^{\mathrm{BW}}=\frac{1}{2}\sqrt{(\mathrm{Re}\{b_{1,n,k}\}-\mathrm{Im}\{a_{1,n,k}\})^{2}+(\mathrm{Im}\{b_{1,n,k}\}+\mathrm{Re}\{a_{1,n,k}\})^{2}}}\\ &{A_{n,k}^{\mathrm{FW}}=\frac{1}{2}\sqrt{(\mathrm{Re}\{b_{1,n,k}\}+\mathrm{Im}\{a_{1,n,k}\})^{2}+(\mathrm{Re}\{a_{1,n,k}\}-\mathrm{Im}\{b_{1,n,k}\})^{2}}}\end{array}
+$$  
+
+and  
+
+$$
+\begin{array}{r l}&{\phi_{n,k}^{0}=\tan^{-1}(\operatorname{Im}\{a_{0,n,k}\}/\operatorname{Re}\{a_{0,n,k}\})}\\ &{\phi_{n,k}^{\mathrm{BW}}=\tan^{-1}[(\operatorname{Im}\{a_{1,n,k}\}-\operatorname{Re}\{b_{1,n,k}\})/(\operatorname{Im}\{b_{1,n,k}\}+\operatorname{Re}\{a_{1,n,k}\})]}\\ &{\phi_{n,k}^{\mathrm{FW}}=\tan^{-1}[(\operatorname{Re}\{b_{1,n,k}\}+\operatorname{Im}\{a_{1,n,k}\})/(\operatorname{Re}\{a_{1,n,k}\}-\operatorname{Im}\{b_{1,n,k}\})]}\end{array}
+$$  
+
+Equation (17) shows that the motion of blade number $i$ for mode number $k$ consists of three components: a symmetric component where all blades deflect simultaneously in blade mode number $n$ with the amplitude $A_{n,k}^{0}$ , and two asymmetric components where the blades deflect with phase shifts of $T/3$ and the amplitudes $A_{n,k}^{\mathrm{BW}}$ and $A_{n,k}^{\mathrm{FW}}$ respectively. As indicated by the superscripts, the asymmetric components represent backward and forward rotor whirling. The direction of the whirl is determined by the sign of the blade-dependent phase shifts $\textstyle{\frac{2}{3}}\pi(i-1)$ .  
+
+The frequency of the symmetric component is $\omega_{k}$ , while the frequencies of the whirling components are shifted by $\pm\Omega$ . The reason for these frequency shifts is that the natural frequency $\omega_{k}$ is given in the ground fixed frame. An observer on the tower top (in the ground fixed frame) will measure the natural frequency $\omega_{k}$ for all three modes. An observer on a blade will also measure the frequency $\omega_{k}$ for symmetric modes, but for backward and forward whirling modes the same observer will measure the frequencies $\omega_{k}+\Omega$ and $\omega_{k}-\Omega$ respectively.  
+
+The damping of the turbine modes is described by the term $\mathrm{e}^{\sigma_{k}t}$ ; the modal damping factor $\sigma_{k}$ determines the exponential decay, or growth, of amplitude for a vibration in mode number $k$ . The computation of the damping factors can therefore be used for stability analyses, e.g. the flutter analysis of a wind turbine.15  
+
+公式(17)表明，对于模态数k，叶片数i的运动由三个分量组成：一个对称分量，所有叶片同时变形，变形模态数为n，振幅为$A_{n,k}^{0}$，以及两个非对称分量，叶片变形具有$T/3$的相位差，振幅分别为$A_{n,k}^{\mathrm{BW}}$和$A_{n,k}^{\mathrm{FW}}$。正如上标所示，非对称分量代表后向和前向风轮摆振。摆振方向由叶片相关的相位差$\textstyle{\frac{2}{3}}\pi(i-1)$的符号决定。
+
+对称分量的频率为$\omega_{k}$，而摆振分量的频率则分别偏移了$\pm\Omega$。产生这些频率偏移的原因是，自然频率$\omega_{k}$是在地面固定参考系中给定的。塔顶观察者（在地面固定参考系中）将测量所有三个模态的自然频率$\omega_{k}$。叶片上的观察者也会测量对称模态的频率$\omega_{k}$，但对于后向和前向摆振模态，相同的观察者将分别测量频率$\omega_{k}+\Omega$和$\omega_{k}-\Omega$ 。
+
+机组模态的阻尼由项$\mathrm{e}^{\sigma_{k}t}$描述；模态阻尼系数$\sigma_{k}$决定了模态数k的振动振幅的指数衰减或增长。因此，阻尼系数的计算可用于稳定性分析，例如风电机组的颤振分析。15
+# Modal Analysis of the Experimental 600 kW Turbine  
+
+This section contains derivations of the natural frequencies and mode shapes of the Bonus $600~\mathrm{kW}$ wind turbine used in the experiments by Thomsen et al.1 A Campbell diagram shows the natural frequencies of the turbine as a function of rotation speed. A mode identification study shows that some of the corresponding mode shapes change dramatically as the rotation speed is changed, owing to modal interactions.  
+
+本节包含对实验中使用的Bonus $600~\mathrm{kW}$ 风电机组的固有频率和模态的推导。一个坎贝尔图（Campbell diagram）显示了风轮的固有频率随转速的变化。**模态识别研究表明，随着转速的变化，一些对应的模态会发生剧烈变化，这是由于模态相互作用的结果**。
+
+# Model Parameters  
+
+The Bonus $600~\mathrm{kW}$ (Mk III C) wind turbine with its $\mathrm{LM}\ 19{\cdot}1\ \mathrm{m}$ blades has a rotor radius of $22\textrm{m}$ with hub and root extension. For this study, two blade modes are chosen $(N=2)$ : first flapwise and first edgewise bending modes for the non-rotating blade. The shape functions $\Phi_{n}^{x}$ and $\Phi_{n}^{y}$ of the entire blade–root–hub assembly are obtained from curve ftis to the results of a finite element model (see Figure 3). The functions are normalized so that the associated modal masses are one. The natural frequencies of the flapwise and edgewise blade modes (as measured in a test stand with the hub and root extension) are 1Ð70 and $2{\cdot}94~\mathrm{Hz}$ respectively.  
+
+The remaining model parameters are obtained by tuning the natural frequencies of the stationary turbine $(\Omega\equiv0)$ ) to the frequencies measured for the particular experimental turbine. The turbine is operating at a speed of about $27.4~\mathrm{rpm}$ $\Omega\approx2{\cdot}87$ rad $\mathrm{s}^{-1}$ ).  
+
+Bonus $600~\mathrm{kW}$ (Mk III C) 风电机组，其 $\mathrm{LM}\ 19{\cdot}1\ \mathrm{m}$ 叶片，风轮半径为 $22\textrm{m}$，包括轴毂和根部延伸。对于本研究，选择了两个叶片模态 $(N=2)$：第一个挥舞模态和第一个摆振模态，用于非旋转叶片。整个叶片–根部–轴毂组件的形状函数 $\Phi_{n}^{x}$ 和 $\Phi_{n}^{y}$，是通过拟合有限元模型的结果获得的（见图3）。这些函数被归一化，使得相关的模态质量为一。在测试支架上测得的挥舞模态和摆振模态的固有频率分别为 1.70 和 $2{\cdot}94~\mathrm{Hz}$。
+
+其余模型参数是通过调整静止机组的固有频率（$\Omega\equiv0$）来匹配特定实验风电机组测得的频率而获得的。风电机组正在以约 $27.4~\mathrm{rpm}$ 的速度运行（$\Omega\approx2{\cdot}87$ rad $\mathrm{s}^{-1}$）。
+![](957376a7e579e755d01141349ef562dc417a67a945cf69141650ed12e2ba5b66.jpg)  
+Figure 3. Mode shape functions for the blade–root–hub assembly: first flapwise mode (left) and first edgewise mode (right). The functions are curve fits to the points obtained from a finite element model with Timoshenko beam elements  
+叶片-根-轮毂组件的模态函数：一阶挥舞模态（左）和一阶摆振模态（右）。这些函数是对有限元模型中Timoshenko梁单元计算所得点的曲线拟合。
+
+Table I. Modelled and measured natural frequencies of the lowest 10 modes of the stationary turbine 
+对固定机组的前10个模态自然频率进行建模和测量。  
+
+
+<html><body><table><tr><td>Mode no.</td><td>Name</td><td>Modelled frequency (Hz)</td><td>Measured frequency (Hz)</td></tr><tr><td>1</td><td>1stshafttorsion</td><td>0.57</td><td>0.56</td></tr><tr><td>2</td><td>1st longitudinal tower bending</td><td>0.73</td><td>0.73</td></tr><tr><td>3</td><td>1st lateral tower bending</td><td>0.78</td><td>0.80</td></tr><tr><td>4</td><td>1st yawing fap</td><td>1·24</td><td>1.20</td></tr><tr><td>5</td><td>1st tilting fap</td><td>1.48</td><td>1.50</td></tr><tr><td>6</td><td>1st symmetric flap</td><td>1.76</td><td>1.76</td></tr><tr><td>7</td><td>1stverticaledgewise</td><td>2.85</td><td>2.90</td></tr><tr><td>8</td><td>1sthorizontaledgewise</td><td>2.95</td><td></td></tr><tr><td>9</td><td>2nd yawing fap</td><td>3.43</td><td>3.60</td></tr><tr><td>10</td><td>2nd tilting flap</td><td>3.74</td><td></td></tr></table></body></html>  
+
+# Natural Frequencies and Mode Shapes  
+
+Modelled and measured natural frequencies of the 10 lowest modes of the stationary turbine are listed in Table I. The measured frequencies are obtained from identification of peaks in frequency responses. For the mode pairs 7/8 and 9/10 it is only possible to identify a single response peak.  
+
+A comparison of modelled and measured frequencies leads to the assumption that the model is sufficiently tuned to the experimental turbine. The naming of the 10 lowest modes is based on animations of their mode shapes and the measured frequency responses. The sequence of modes is typical for a turbine of this size, with the exception that the 1st shaft torsion may be higher than the tower bending modes. The longitudinal tower bending always lies slightly lower than the lateral tower bending mode, because it contains some tilting of the rotor, which has a large inertia.  
+
+Modes 4 and 5 are the 1st tilt/yaw modes involving the flapwise blade mode. The yaw mode most often lies lower than the tilt mode, because towers are stiffer in tilt than in yaw. The sixth mode is the 1st symmetrical flap mode, where the blades vibrate simultaneously in the flapwise blade mode in counter-phase with a longitudinal tower vibration. This coupling causes the symmetric flap frequency to lie slightly above the frequency of the flapwise blade mode.  
+
+Modes 7 and 8 involve the edgewise blade mode, and their frequencies are close to the edgewise blade frequency of $2{\cdot}94~\mathrm{Hz}$ (with an average slightly below). In both modes the blades vibrate edgewise against each other so that they cancel out the torsional moment at the rotor centre. The two modes differ in the direction of the reactive force at the rotor centre, as indicated by their names: 1st vertical and 1st horizontal edgewise mode. The sequence of these two modes is given by the vertical and horizontal stiffness of the rotor support.  
+
+The last two modes listed in Table I are the 2nd tilt/yaw modes, where the rotor blades are tilting and yawing in counter-phase with the tilt and yaw of the nacelle. Again the yaw mode lies below the tilt mode because of the lower yaw than tilt stiffness of the tower.  
+
+Figure 4 shows how the natural frequencies $\omega_{k}$ of these 10 lowest modes change with the rotation speed $\Omega$ from standstill to the operation speed. The natural frequencies of the tower bending modes and the shaft torsion mode are constant with rotation speed. The frequency of the symmetric flap mode increases owing to centrifugal stiffening of flapwise bending. The natural frequencies of the asymmetric rotor modes change with rotation speed owing to gyroscopic effects.  
+
+The naming of the modes at operation listed in the table in Figure 4 requires some introduction. All pairs of asymmetric rotor modes at standstill (modes 4/5, 7/8 and 9/10) become pairs of rotor whirling modes owing to the rotation, e.g. the 1st tilt/yaw modes become the 1st flapwise whirling modes. The frequencies of the backward whirling (BW) modes decrease with rotation speed, whereas the frequencies of the forward whirling (FW) modes increase with rotation speed.  
+
+表I列出了静止机组前10阶模态的固有频率，这些频率是通过识别频率响应中的峰值获得的。对于7/8阶和9/10阶模态对，仅能识别出一个响应峰值。
+
+对模拟和测量频率进行比较，表明模型已经足够地针对实验机组进行了调整。前10阶模态的命名是基于其模态形状的动画和测量频率响应。这些模态的顺序对于这种尺寸的风电机组来说是典型的，但1阶主轴扭转可能高于塔架弯曲模态。纵向塔架弯曲总是略低于横向塔架弯曲模态，因为横向塔架弯曲包含一些叶片摆动，而叶片具有很大的惯性。
+
+第4阶和第5阶模态是第1阶倾转/偏航模态，涉及挥舞叶片模态。偏航模态通常低于倾转模态，因为塔架在倾转方向上比在偏航方向上更刚。第6阶模态是第1阶对称挥舞模态，其中叶片以逆相位与纵向塔架振动同时在挥舞方向上振动。这种耦合导致对称挥舞频率略高于挥舞叶片模态的频率。
+
+第7阶和第8阶模态涉及摆振叶片模态，它们的频率接近于摆振叶片频率的$2{\cdot}94~\mathrm{Hz}$（平均值略低于此值）。在这两个模态中，叶片相互作用，在摆振方向上振动，从而抵消了风轮中心的扭转力矩。这两个模态的区别在于风轮中心反作用力的方向，这由它们的名称所示：第1阶垂直摆振模态和第1阶水平摆振模态。这两个模态的顺序由风轮支撑的垂直和水平刚度决定。
+
+表I中列出的最后两个模态是第2阶倾转/偏航模态，其中风轮叶片以与机舱倾转和偏航相反的相位倾转和偏航。同样，由于塔架在偏航方向上的刚度低于倾转方向，偏航模态低于倾转模态。
+
+图4显示了这些前10阶简正模态的固有频率$\omega_{k}$随转速$\Omega$从静止到额定转速的变化情况。塔架弯曲模态和主轴扭转模态的固有频率随转速变化而保持不变。由于挥舞方向的离心加固，对称挥舞模态的固有频率增加。由于陀螺效应，非对称风轮模态的固有频率随转速变化。
+
+图4中列出的额定转速下模态的命名需要一些介绍。所有在静止状态下的非对称风轮模态对（模态4/5、7/8和9/10）都因转动而变成风轮回转模态对，例如，第1阶倾转/偏航模态变成第1阶挥舞回转模态。向后回转（BW）模态的频率随转速降低，而向前回转（FW）模态的频率随转速升高。
+
+
+![](1668404ad97d76291f188361e7728ca1d61b3c36e9a82d524b6e3b4ca4ff9ea2.jpg)  
+Figure 4. Campbell diagram: natural frequencies as a function of rotation speed of the rotor. The points denote the measured natural frequencies for the turbine at standstill  
+
+This splitting of the natural frequencies of a whirling mode pair is related to the co-ordinate system of observation through gyroscopic effects. The natural frequencies $\omega_{k}$ in Figure 4 are observed from the ground fixed frame. Equation (17) shows that in a co-rotating blade frame the frequency of a symmetric rotor mode remains $\omega_{k}$ , whereas the frequencies of BW and FW modes become $\omega_{k}+\Omega$ and $\omega_{k}-\Omega$ respectively. If only blade mode number $n$ is involved in the whirling modes, their frequencies split in the ground fixed frame about the natural frequency $\omega_{n}$ of this blade mode. In this ideal case the observer on the blade will measure the same frequency $\omega_{n}=\omega_{k}^{\mathrm{BW}}+\Omega=\omega_{k}^{\mathrm{FW}}-\Omega$ for both the BW and FW modes.  
+
+This ideal condition is affected by structural asymmetry of the turbine and coupling of blade modes in the turbine modes. Such modal interactions may occur when the natural frequencies of two turbine modes come close. Figure 4 shows that the frequencies of the 1st FW edgewise mode (mode 8) and the 2nd BW flapwise mode (mode 9) become close at about $\Omega=2{\cdot}5$ rad $\mathrm{s}^{-1}$ . These two modes interact, which can be shown by the flapwise and edgewise whirling components in their mode shapes.  
+
+Figure 5 shows the FW and BW components of the flapwise and edgewise blade modes in modes 7–10. These whirling components are computed from the eigenvectors in multi-blade co-ordinates using equation (18). The dominating amplitudes for modes 7 and 10 show that these modes are respectively a BW edgewise and an FW flapwise mode. However, there are no dominant modal amplitude for modes 8 and 9 over the whole range of rotation speeds. It seems that these modes interchange mode shapes: mode 8 can be defined as an FW edgewise mode and mode 9 as a BW flapwise mode at rotation speeds below 2Ð5 rad $\mathrm{s}^{-1}$ , and vice versa above this speed.  
+
+At the operation speed, mode 9 must be characterized as the 1st FW edgewise mode with some content of flapwise whirling. This content is now suggested to explain the observed difference in aerodynamic damping of the BW and FW edgewise modes (modes 7 and 9).  
+
+# Why the Measured Difference in Aerodynamic Damping?  
+
+Thomsen et al.1 have estimated the total damping of the two edgewise whirling modes (see Figure 1). The results show that the FW edgewise mode (mode 9) is more damped than the BW edgewise mode (mode 7). The structural damping of these two modes is assumed to be the same, because their natural frequencies and mode shapes are almost identical. Thus the difference in total damping is assumed to be caused by a difference in aerodynamic damping.  
+
+![](b27d093788eedd61e8374571ef2a21361ad23a35d29ba894deefacc295c3a9c7.jpg)  
+Figure 5. Flapwise and edgewise whirling components of the modal blade amplitudes for modes 7–10 computed from the eigenvectors using equation (18)  
+
+The aim of the experiment was to measure the damping of edgewise blade vibrations during turbine operation by exciting the corresponding turbine modes using an exciter mounted in the nacelle. The presented theoretical modal analysis of the turbine was not available for the experiment. The excitation of the edgewise whirling modes was therefore based on tuning the excitation frequency to obtain the maximum edgewise blade response. The two initial guesses, however, were based on the assumption that the ideal condition for frequency splitting existed, i.e. the excitation frequency should be the natural frequency of the edgewise blade mode plus or minus the operation speed (1P).  
+
+Figure 6 shows a zoom in the Campbell diagram (Figure 4) on the natural frequencies of modes 7–9, together with lines of the $\pm1\mathrm{P}$ splitting about the edgewise blade frequency $(2{\cdot}94\ \mathrm{Hz})$ . These lines intersect with the natural frequencies of modes 7 and 9 at the operation speed, showing that the initial guesses on the experimental excitation frequencies are close. However, the modal analysis has also shown that mode 9 is not a pure edgewise whirling mode; it has some flapwise component. This component affects the direction in which the blades are vibrating relative to the rotor plane.  
+
+# Effective Direction of Blade Vibration  
+
+Figure 7 shows how the blade cross-section at $90\%$ radius is moving in and out of the rotor plane during vibrations in the two edgewise whirling modes (modes 7 and 9). The motion of the cross-section includes the motion of the rotor support, i.e. the shaft and tower deformations add to the effective blade vibration. The traces show that the blades move more out of the rotor plane in the FW edgewise mode than in the  
+
+![](fbc1b019314328084a7e66eb9b851b544a45506606be095b8746d98395d8f071.jpg)  
+Figure 6. A zoom in the Campbell diagram (Figure 4) on the natural frequencies of modes 7–9. The broken lines show the $\pm\Omega$ splitting about the edgewise blade frequency $\left(2{\cdot}94\,H z\right)$  
+
+![](b8ac1f69e8e67bd89e848dd13bf3bc4bab5b265f8838cf15dc851116eb006242.jpg)  
+Figure 7. Motion of a blade cross-section at $90\%$ radius for the 1st BW edgewise mode (left) and the 1st FW edgewise mode (right). The motion includes the tower and shaft deformations  
+
+BW edgewise mode. This difference in blade vibration for the two modes can explain why Thomsen et al.   
+observed that the FW mode is more damped than the BW mode.  
+
+A quasi-steady aerodynamic analysis in the Appendix shows that the aerodynamic damping of a blade cross-section is lowest for vibrations close to the rotor plane. It is presumed in the analysis that the blade cross-section is performing a small elliptical motion, with the major axes being parallel and perpendicular to the rotor plane. This type of motion is assumed to characterize the qualitative pattern of the traces in Figure 7, except for a slight tilt of the major axes which originates from the direction of blade vibration for the edgewise blade mode (see Figure 3). An effective direction of blade vibration is defined for the elliptical motion of a blade cross-section as  
+
+$$
+\theta_{\mathrm{eff}}=\tan^{-1}\left({\frac{\mathrm{max.~amplitude~out~of~rotor~plane}}{\mathrm{max.~amplitude~in~rotor~plane}}}\right)
+$$  
+
+Because the elliptical motion is assumed to have the major axes parallel and perpendicular to the rotor plane, the maximum amplitudes are positive, yielding $0^{\circ}<\theta_{\mathrm{eff}}<90^{\circ}$ . Using this definition on the motion of the cross-section at $90\%$ radius in Figure 7, it is found that $\theta_{\mathrm{eff}}\approx9^{\circ}$ and $30^{\circ}$ for the 1st BW and FW edgewise modes respectively. Computations of $\theta_{\mathrm{eff}}$ along the blade show that the entire blade is vibrating more out of the rotor plane in the 1st FW edgewise mode than in the 1st BW edgewise mode. Hence, based on quasisteady aerodynamics, this behaviour can explain the measured difference in aerodynamic damping of these two modes.  
+
+# Improved Modal Dynamics  
+
+The out-of-plane motion of the blades in the 1st FW edgewise mode is mainly due to a component of the flapwise blade mode through the previously described modal interaction with the 2nd BW flapwise mode (see Figure 5). The 1st BW edgewise mode is not interacting with a flapwise whirling mode; the blades are mainly vibrating in the edgewise blade mode (with the additional motion due to the tower and shaft deformations). Is it possible to add a flapwise component to the BW edgewise mode without removing the flapwise component in the FW edgewise mode, thereby increasing the overall aerodynamic damping of the turbine?  
+
+An example shows that such improvement of the modal dynamics of the Bonus $600~\mathrm{kW}$ turbine may be possible. The natural frequency of the edgewise blade mode is chosen to be the design parameter in this optimization problem. Note that in reality an improved design must be the result of an integrated process towards a ‘global optimum’, where not only a single design parameter (the edgewise blade stiffness) is considered to achieve a single objective (the effective direction of blade vibration).  
+
+Figure 8 shows the effective directions of blade vibration defined by equation (20) at the $90\%$ radius for modes 7–10 as function of the edgewise blade frequency. The modes are here named by their numbers, because their modal interactions prevent a unique mode identification over the entire range of design variation. These modes are all interesting, because they involve the edgewise blade mode. Shaft torsion and lateral tower bending modes also involve blade vibrations close to the rotor plane; however, the resulting low aerodynamic damping is compensated by high structural damping due to the generator slip of an asynchronous machine.  
+
+A local optimum is marked at an increase of $22\%$ in the edgewise blade frequency, from 2Ð94 to $3{\cdot}59\ \mathrm{Hz}$ . For this design the effective direction of blade vibration for mode 7 has tripled to approximately $27^{\circ}$ . The effective direction of blade vibration has also increased for mode 8; however, the blades are now vibrating closer to the rotor plane in modes 9 and 10. In fact, mode 10 would be the most critical mode for higher edgewise blade frequencies.  
+
+![](e36f60f30c18c521b819097fdc94747ac3d39dfc9c912bc63f4000e25e395b40.jpg)  
+Figure 8. Effective directions of blade vibration at $90\%$ radius with respect to the rotor plane for modes 7–10 as a function of the edgewise blade frequency  
+
+The improved design (with respect to the effective direction of blade vibration) uses the modal interaction with 2nd flapwise whirling modes to increase the flapwise component in the 1st edgewise whirling modes. It is not a coincidence that the optimal edgewise blade frequency is almost identical to the average frequency of the 2nd tilt and yaw modes at standstill (see Table I and Figure 4). Figure 9 shows a Campbell diagram for modes 7–10 of the improved turbine, together with the $\pm1P$ splitting about the edgewise blade frequency. The ideal condition with a single blade mode involved in the rotor whirling mode is not present, and the modal interaction of edgewise and flapwise whirling is complete.  
+
+For this particular turbine it may be impossible to increase the edgewise blade frequency by $22\%$ ; however, by changing other design parameters, it may be possible to lower the frequencies of the 2nd flapwise whirling modes to obtain the desired modal interaction. A rule of thumb in structural design for maximum out-ofrotor-plane blade vibration can be formulated as: Make sure the 1st edgewise blade frequency lies between the frequencies of the 2nd tilt and yaw modes. Whether and how this design guideline can be followed will depend on the wind turbine.  
+
+# Conclusion  
+
+A theoretical modal analysis of a three-bladed wind turbine shows that the blades cannot be considered as separate dynamic components of the turbine. Blade vibrations in the different modes of the turbine are strongly affected by the dynamics of the shaft, nacelle and tower. It is therefore necessary to consider the dynamics of the complete turbine when estimating the risk of stall-induced blade vibrations.  
+
+The analysis of a stall-controlled $600~\mathrm{kW}$ turbine suggests why the forward and backward edgewise whirling modes have been observed to have different aerodynamic damping.1 The explanation is a difference in the blade vibration for the two edgewise whirling modes. The forward whirling mode interacts with a flapwise whirling mode, whereby the blades vibrate more out of the rotor plane, yielding an increased aerodynamic damping.  
+
+Variations of the turbine design indicate that the structural dynamics of a wind turbine can be tailored to obtain a higher aerodynamic damping of its critical modes. A rule of thumb in structural design to lower the risk of stall-induced edgewise vibrations can be formulated as: Make sure the 1st edgewise blade frequency lies between the frequencies of the 2nd tilt and yaw modes. Whether and how this design guideline can be followed will depend on the wind turbine.  
+
+![](43f05a925486cdde7941cb8be05a93435038eb584261ba1bd5fff97bffbf8dc4.jpg)  
+Figure 9. Natural frequencies of modes 7–10 for the improved turbine. The broken lines show the $\pm\Omega$ splitting about the new edgewise blade frequency of $3{\cdot}59\,H z$  
+
+# Acknowledgements  
+
+This work was partly funded by the Danish Energy Agency under ENS 1363/00-0006, which is gratefully acknowledged. The author would also like to thank Kenneth Thomsen and Jørgen Thirstrup Petersen for their many invaluable discussions and suggestions.  
+
+# Appendix: Aerodynamic Damping  
+
+This appendix deals with the aerodynamic damping of blade vibrations based on quasi-steady aerodynamics. Unsteady aerodynamic effects can be dominating when the blades operate in stall; however, the basic qualitative behaviour of the aerodynamic damping is assumed to be described by quasi-steady aerodynamics (see e.g. References 8 and 16 for details on dynamic stall effects).  
+
+The aim is to determine the aerodynamic damping of a vibrating blade cross-section. The cross-section is assumed to translate in and out of the rotor plane with the velocities $\dot{u}_{x}$ and $\dot{u}_{y}$ respectively, as shown in Figure 10. Torsional motion is neglected. The velocity triangle shows that the relative flow angle with respect to the rotor plane and the relative flow speed are  
+
+$$
+\phi=\tan^{-1}\left(\frac{V-\dot{u}_{y}}{r\Omega+\dot{u}_{x}}\right),\qquad W=\sqrt{(V-\dot{u}_{y})^{2}+(r\Omega+\dot{u}_{x})^{2}}
+$$  
+
+where $V$ is the wind speed, and the product of radius and rotation speed, $r\Omega$ , is the steady tangential speed of the cross-section. Note that induced velocities from the wake behind the turbine are neglected, which has no qualitative influence on the derivations. The angle of attack is the difference between the relative flow angle and the twist of the cross-section: $\alpha=\phi-\theta$ .  
+
+The assumption of quasi-steady conditions (neglecting added mass effects) means that the lift and drag on the cross-section can be written as  
+
+$$
+L(\dot{u}_{x},\dot{u}_{y})=\frac{1}{2}c\rho W^{2}C_{\mathrm{L}}(\alpha),\quad D(\dot{u}_{x},\dot{u}_{y})=\frac{1}{2}c\rho W^{2}C_{\mathrm{D}}(\alpha)
+$$  
+
+where $c$ is the chord length, $\rho$ is the air density, and $C_{\mathrm{{L}}}$ and $C_{\mathrm{D}}$ are the steady aerodynamic lift and drag coefficients given as functions of the angle of attack. The lift and drag forces are non-linear functions of the velocities $\dot{u}_{x}$ and $\dot{u}_{y}$ of the cross-section.  
+
+Assuming that these velocities are small compared with the steady flow velocities $\dot{u}_{x}<<r\Omega$ and $\dot{u}_{y}<<V$ ), the aerodynamic forces (22) are linearized and projected onto the co-ordinate system of the rotor plane $(x,\,y)$ . The resulting linear quasi-steady aerodynamic forces can be written as  
+
+$$
+\mathrm{F}=\left\{\!\!\begin{array}{c}{F_{x}}\\ {F_{y}}\end{array}\!\right\}=-\frac{1}{2}c\rho W_{0}\left[\!\!\begin{array}{c c}{c_{x x}}&{c_{x y}}\\ {c_{y x}}&{c_{y y}}\end{array}\!\!\right]\left\{\!\!\begin{array}{c}{{\dot{u}}_{x}}\\ {{\dot{u}}_{y}}\end{array}\!\!\right\}
+$$  
+
+![](44d5a13a9fc7e28b57431322ff87fd24734d5f686e5cdd783fa2c62cb397a129.jpg)  
+Figure 10. Velocity triangle for a blade cross-section (without induced velocities from the turbine wake)  
+
+where $W_{0}=\sqrt{V^{2}+r^{2}\Omega^{2}}$ is the steady relative velocity, and the diagonal elements of the damping coefficient matrix are  
+
+$$
+\begin{array}{r l}&{c_{x x}=(3C_{\mathrm{D}}+C_{\mathrm{L}}^{\prime}+(C_{\mathrm{D}}-C_{\mathrm{L}}^{\prime})\cos2\phi_{0}-(C_{\mathrm{L}}+C_{\mathrm{D}}^{\prime})\sin2\phi_{0})/2}\\ &{c_{y y}=(3C_{\mathrm{D}}+C_{\mathrm{L}}^{\prime}-(C_{\mathrm{D}}-C_{\mathrm{L}}^{\prime})\cos2\phi_{0}+(C_{\mathrm{L}}+C_{\mathrm{D}}^{\prime})\sin2\phi_{0})/2}\end{array}
+$$  
+
+where $()^{\prime}\equiv\mathrm{d}/\mathrm{d}\alpha$ , and $\phi_{0}=\tan^{-1}(V/r\Omega)$ is the steady flow angle. The lift and drag coefficients and their derivatives are evaluated at the steady angle of attack $\alpha_{0}=\phi_{0}-\theta$ .  
+
+Previous studies of quasi-steady aerodynamic damping have presented linear aerodynamic forces identical to equation (23), although their expressions for the damping coefficient have another form (see e.g. References 5, 7 and 8). These studies deal with the damping of blade vibrations, where the blade cross-sections move in straight lines in the $(x,\,y)$ plane. This is the case for a blade mounted in a test stand; however, the present modal analysis of the rotating turbine shows that cross-sections move in a more elliptical motion when the blades are vibrating in a turbine mode (see Figure 7).  
+
+It is assumed that the cross-section is moving in an ellipse with the major axes parallel and perpendicular to the rotor plane, as shown in Figure 11. Using definition (20) of the effective direction of vibration $\theta_{\mathrm{eff}}$ to describe the ratio of the major axes, this elliptical motion of a blade cross-section can be described by the vector  
+
+$$
+\mathbf{u}=a\left\{\begin{array}{c}{\cos{\theta_{\mathrm{eff}}}\cos{\omega t}}\\ {\sin{\theta_{\mathrm{eff}}}\sin{\omega t}}\end{array}\right\}
+$$  
+
+where $a$ and $\omega$ are the amplitude and frequency of the elliptical motion.  
+
+The aerodynamic damping can be approximated as the negative work done by the aerodynamic forces on the blade cross-section over one period of oscillation, $T\equiv2\pi/\omega$ . Positive aerodynamic work corresponds to negative aerodynamic damping. For the elliptical motion (25) the aerodynamic work can be derived as  
+
+$$
+W_{\mathrm{aero}}=\int_{0}^{T}\mathrm{F^{T}i d}t=-\frac{\pi}{2}\omega a^{2}c\rho W_{0}c_{\mathrm{eff}}
+$$  
+
+where $c_{\mathrm{eff}}$ is an effective aerodynamic damping coefficient given by  
+
+$$
+c_{\mathrm{eff}}=\frac{1}{2}(c_{x x}+c_{y y})+\frac{1}{2}(c_{x x}-c_{y y})\cos{2\theta_{\mathrm{eff}}}
+$$  
+
+When this coefficient is negative, the aerodynamic damping of the elliptical motion is negative. Note that the elliptical motion (25) is clockwise; the work of the linear aerodynamic forces is identical to (26) for counter-clockwise motion (this is not the case when unsteady aerodynamic effects are included).  
+
+Expression (26) shows that $c_{\mathrm{eff}}=c_{x x}$ for purely in-plane motion of the cross-section $(\theta_{\mathrm{eff}}=0^{\circ}$ ), whereas $c_{\mathrm{eff}}=c_{\mathrm{yy}}$ for purely out-of-plane motion of the cross-section $(\theta_{\mathrm{eff}}=90^{\circ}$ ). For elliptical motion between purely in- and out-of-plane motion the effective damping coefficient lies between $c_{x x}$ and $c_{y y}$ . Figure 12 shows these two damping coefficients for the cross-section at $90\%$ radius as a function of operational wind speeds of the experimental turbine analysed in this article. The grey area shows the effective damping of elliptical motion of the cross-section. Note that $c_{x x}<c_{y y}$ at all operational wind speeds, which shows that the aerodynamic damping at this radius is always lowest for the blade vibration closest to the rotor plane.  
+
+![](3175162a5b5d2b2ebe20ea696e4b39ed72755d768268afb5258e86f9e462fa61.jpg)  
+Figure 11. Elliptic motion of the blade cross-section described by the vector u (equation (25)). The effective direction of blade vibration is given by $\theta_{\mathrm{eff}}$  
+
+![](0b3ecfb643a49c2f849e91f75191baec0437423ee2f87cd11b7b34edcffe9c73.jpg)  
+Figure 12. Quasi-steady aerodynamic damping coefficients (equation (24)) for in-plane and out-of-plane vibrations of the blade cross-section at $90\%$ radius for the experimental $600\,k W$ turbine.1 Effective damping coefficients for elliptical motion lie in the grey area between the curves  
+
+The simple expression (26) enables a qualitative understanding of the critical parameters for stall-induced vibrations. The mean $(c_{x x}+c_{y y})/2$ determines a mean aerodynamic damping of arbitrary blade vibrations. The variation $(c_{x x}-c_{y y})/2$ about this mean determines how the damping varies with the effective direction of vibration in and out of the rotor plane. Using equation (24), the sum and difference of the two aerodynamic coefficients become  
+
+$$
+\begin{array}{r l}&{c_{x x}+c_{y y}=3C_{\mathrm{D}}+C_{\mathrm{L}}^{\prime}}\\ &{c_{x x}-c_{y y}=(C_{\mathrm{D}}-C_{\mathrm{L}}^{\prime})\cos2\phi_{0}-(C_{\mathrm{L}}+C_{\mathrm{D}}^{\prime})\sin2\phi_{0}}\end{array}
+$$  
+
+which shows that the drag, and the lift gradient add to the mean aerodynamic damping, e.g. a negative lift gradient decreases the mean damping. It also shows that the variation of damping with the direction of vibration depends on all aerodynamic profile coefficients weighted with the steady relative flow angle $\phi_{0}$ . It can be deduced from (28) that the aerodynamic damping of in-plane motion increases proportionally to drag but decreases proportionally to lift, lift gradient and drag gradient, because $0^{\circ}\le\phi_{0}\le90^{\circ}$ for most operation conditions.  
+
+# References  
+
+1. Thomsen K, Petersen JT, Nim E, Øye S, Petersen B. A method for determination of damping for edgewise blade vibrations. Journal of Wind Energy 2001; 3: 233–246.   
+2. Hansen MH, Thomsen K, Petersen JT. Rotor whirling modes and the relation to their aerodynamic damping. 2001 European Wind Energy Conference, Copenhagen, 2001; 422–425.   
+3. Rasmussen F, Petersen JT, Winkelaar D, Rawlinson-Smith R. Response of stall regulated wind turbines—stall induced vibrations. Technical Report Risø-R-691(EN), Risø National Laboratory, Roskilde, 1993.   
+4. Stiesdal H. Extreme wind loads on stall regulated wind turbines. In Proceedings of the 16th British Wind Energy Association Conference, Elliot G (ed.). Mechanical Engineering Publications Limited: London; 1994, 101–106.   
+5. Bj¨orck A, Dahlberg J, ¨Ostman A, Ganander H. Computations of aerodynamic damping for blade vibrations in stall.   
+1997 European Wind Energy Conference, Dublin, 1997; 503–507.   
+6. Petersen JT, Thomsen K, Madsen HAa. Local blade whirl and global rotor whirl interaction. Technical Report Risø- R-1067(EN), Risø National Laboratory, Roskilde, 1998.   
+7. Petersen JT, Madsen HAa, Bj¨orck A, Enevoldsen P, Øye S, Ganander H, Winkelaar D. Prediction of dynamic loads and induced vibrations in stall. Technical Report Risø-R-1045(EN), Risø, Roskilde, 1998.   
+8. Rasmussen F, Petersen JT, Madsen HAa. Dynamic stall and aerodynamic damping. ASME Journal of Solar Energy Engineering 1999; 121: 150–155.   
+9. Petersen JT, Thomsen K, Madsen HAa. Stall strips can control edgewise vibrations. Technical Report AED–RB–6 (EN), Risø National Laboratory, Roskilde, 1998.   
+10. Anderson C, Heerkes H, Yemm R. The use of blade-mounted dampers to eliminate edgewise stall vibration. 1999 European Wind Energy Conference, Nice, 1999; 207–211.   
+11. Petersen JT. The aeroelastic code HAWC—model and comparisons. In State of the Art of Aeroelastic Codes for Wind Turbine Calculations, vol. Annex XI, Pedersen BM (ed). International Energy Agency/Technical University of Denmark: Lyngby, 1996; 129–135.   
+12. Johnson W. Helicopter Theory. Princeton University Press: Princeton, NJ, 1980.   
+13. Peters DA. Fast Floquet theory and trim for multi-bladed rotorcraft. Journal of the American Helicopter Society 1994;   
+39(4): 82–89.   
+14. Ziegler H. Principles of Structural Stability (2nd edn). Basel and Stuttgart, Birkh¨auser, 1977.   
+15. Hansen MH. Vibrations of a three-bladed wind turbine due to classical flutter. 2002 ASME Wind Energy Symposium, Reno, NV, 2002.   
+16. Leishman JG. Principles of Helicopter Aerodynamics. Cambridge University Press: Cambridge, 2000.