diff --git a/书籍/力学书籍/CASEstab_theory_manual/auto/CASEstab_theory_manual.md b/书籍/力学书籍/CASEstab_theory_manual/auto/CASEstab_theory_manual.md new file mode 100644 index 0000000..e6726d3 --- /dev/null +++ b/书籍/力学书籍/CASEstab_theory_manual/auto/CASEstab_theory_manual.md @@ -0,0 +1,1239 @@ +# Theory manual for CASEStab + +Morten Hartvig Hansen + +A computational-aero-servo +elastic simulation and +stability tool for wind +turbines + +# Colophon + +Theory manual for the hydro­aero­servo­elastic software package CASEStab This report ... + +# Contents + +1 Structural equations of motion 4 + +1.1 Lagrange equations of articulated substructures 4 +1.1.1 Generic topology of structure . 5 +1.1.2 Inertia forces on and from a substructure 7 +1.1.3 Gravitational forces on and from a substructure . 12 +1.1.4 Linear elastic forces inside a substructure . 12 +1.1.5 Linear purely dissipative forces inside a substructure 12 +1.1.6 Generalized forces from external forces and moments 12 +1.2 Substructure models . 12 +1.2.1 Rigid body substructure . 12 +1.2.2 Co­rotational beam substructure 12 +1.2.3 Linear super­element substructure 16 + +# 2 Rotor aerodynamics 17 + +2.1 Generalized aerodynamic forces on and from a blade 17 +2.1.1 Generalized aerodynamic forces on and from a substructure of the blade . 19 +2.2 Aerodynamic forces and moment at a cross­section . 20 +2.2.1 Velocity triangles . . . . 21 +2.2.2 Unsteady aerodynamic model 22 +2.3 Aerodynamic discretization of a blade 23 +2.4 Wake induction model . . 24 +2.4.1 Stationary steady­state rotor aerodynamics . . 24 +2.4.2 Periodic steady­state rotor aerodynamics . 24 + +# 3 Aeroelastic equations of motion 25 + +3.1 Generalized aerodynamic forces for different substructure types . 25 +3.1.1 Rigid body substructure . . . . . 26 +3.1.2 Co­rotational beam substructure . . 26 + +# A Co­rotational beam element method 33 + +A.1 Nonlinear kinematic formulation . 33 +A.1.1 Element coordinate system 33 +A.1.2 Elastic stiffness matrix . 33 +A.2 Model input for co­rotational beam structures 33 +A.3 Structural HAWC2 input translator 33 +A.4 Cross­sectional transformation 34 + +B Aerodynamic data input and processing 36 + +B.1 Input files 36 + +C Dynamic stall model data 37 + +D Wind field 38 + +E Mathematics 39 + +E.1 Definitions . 39 +E.2 Finite Rotations . 39 +E.3 Miscellaneous formulas 39 + +# 1 Structural equations of motion + +This chapter contains the derivations of the nonlinear equations of motion based on the kinematic formulation of articulated substructures to form the entire turbine structure. + +# 1.1 Lagrange equations of articulated substructures + +The Lagrangian of a structure ${\cal{L}}=T-V$ is given by its total kinetic energy $T$ and the total potential energy $V$ of the conservative forces acting on the structure, e.g. gravity and elastic forces. Using Lagrange’s equations, the nonlinear equations of motion can be derived as + +$$ +{\frac{d}{d t}}\left({\frac{\partial L}{\partial{\dot{q}}_{i}}}\right)-{\frac{\partial L}{\partial q_{i}}}+{\frac{\partial D}{\partial{\dot{q}}_{i}}}=Q_{i}\;\;{\mathsf{f o r}}\;\;i=1,\ldots,N_{D} +$$ + +where $D$ is the Rayleigh dissipation function used to model the internal energy dissipation, and $Q_{i}$ is the gen­ eralized force for the generalized coordinate $q_{i}$ due to non­conservative forces handled in the next sections. + +The Lagrangian is a function of the displacements, velocities, and time + +$$ +L=T\left(t,\mathbf{q},{\dot{\mathbf{q}}}\right)-V\left(t,\mathbf{q}\right) +$$ + +where the potential energy of the conservative forces are independent of velocities. Using (1.2) and the chain rule, the first term of (1.1) can be expanded as + +$$ +\sum_{j=1}^{N_{D}}\frac{\partial^{2}T}{\partial{\dot{q}}_{i}\partial{\dot{q}}_{j}}{\ddot{q}}_{j}+\sum_{j=1}^{N_{D}}\frac{\partial^{2}T}{\partial{\dot{q}}_{i}\partial q_{j}}{\dot{q}}_{j}+\frac{\partial}{\partial t}\bigg(\frac{\partial T}{\partial{\dot{q}}_{i}}\bigg)-\frac{\partial T}{\partial q_{i}}+\frac{\partial D}{\partial{\dot{q}}_{i}}+\frac{\partial V}{\partial q_{i}}=Q_{i} +$$ + +for $i=1,2,\dots,N_{D}$ . Here, the first term constitutes the acceleration dependent forces, whereas the other terms are only dependent on the time $t$ and the state­variables, the displacements $\mathbf{q}$ and velocities $\dot{\mathbf{q}}$ . We will now take a closer look at the first four inertia force terms given by the kinetic energy. + +The total kinetic energy is the integral of the kinetic energy of each particle over the entire volume $\mathcal{V}$ of the structure: + +$$ +T=\int_{\mathcal{V}}\frac{1}{2}\;\rho\;\dot{\mathbf{r}}^{T}\dot{\mathbf{r}}\;d\mathcal{V} +$$ + +where $()^{T}$ denotes to the transpose of a matrix or a vector (single columned matrix), and r˙ is the velocity vector of the particle given as the time derivative of its position vector $\boldsymbol{\mathsf{r}}=\boldsymbol{\mathsf{r}}(t,\mathbf{q})$ that may be explicit time­dependent e.g. for sub­structures that are rotating with a prescribed average speed. The velocity vector can be expanded to + +$$ +\dot{\pmb{r}}=\frac{d\pmb{r}}{d t}=\sum_{j=1}^{N_{D}}\frac{\partial\pmb{r}}{\partial q_{j}}\dot{q}_{j}+\frac{\partial\pmb{r}}{\partial t} +$$ + +from which these properties of the position and velocity vectors can be shown + +$$ +{\frac{\partial{\dot{\mathbf{r}}}}{\partial{\dot{q}}_{i}}}={\frac{\partial\mathbf{r}}{\partial q_{i}}}\quad{\mathsf{a n d}}\quad{\frac{\partial^{2}{\dot{\mathbf{r}}}}{\partial{\dot{q}}_{i}\partial q_{j}}}={\frac{\partial\mathbf{r}}{\partial q_{i}\partial q_{j}}}={\frac{\partial^{2}{\dot{\mathbf{r}}}}{\partial q_{i}\partial{\dot{q}}_{j}}} +$$ + +Substituting (1.4) into (1.3) and using the properties of (1.5) and (1.6), the nonlinear equations of motion can be written as + +$$ +\sum_{j=1}^{N_{D}}m_{i j}\ddot{q}_{j}+\sum_{j=1}^{N_{D}}g_{i j}\dot{q}_{j}+F_{c,i}+\sum_{j=1}^{N_{D}}\sum_{k=1}^{N_{D}}h_{i j k}\dot{q}_{j}\dot{q}_{k}+\frac{\partial V}{\partial q_{i}}+\frac{\partial D}{\partial\dot{q}_{i}}=Q_{i} +$$ + +where the summation coefficients and the acceleration force $F_{c,i}$ can be derived solely from the partial time and displacement derivatives of the position vector as + +$$ +\begin{array}{l}{\displaystyle m_{i j}=\int_{\mathcal{V}}\rho\left(\frac{\partial\mathbf{r}^{T}}{\partial q_{i}}\frac{\partial\mathbf{r}}{\partial q_{j}}\right)d\boldsymbol{\nu}}\\ {\displaystyle g_{i j}=\int_{\mathcal{V}}2\rho\left(\frac{\partial\mathbf{r}^{T}}{\partial q_{i}}\frac{\partial^{2}\mathbf{r}}{\partial t\partial q_{j}}\right)d\boldsymbol{\nu}}\\ {\displaystyle h_{i j k}=\int_{\mathcal{V}}\rho\left(\frac{\partial\mathbf{r}^{T}}{\partial q_{i}}\frac{\partial^{2}\mathbf{r}}{\partial q_{j}\partial q_{k}}\right)d\boldsymbol{\nu}}\\ {\displaystyle F_{c,i}=\int_{\mathcal{V}}\rho\left(\frac{\partial\mathbf{r}^{T}}{\partial q_{i}}\frac{\partial^{2}\mathbf{r}}{\partial t^{2}}\right)d\boldsymbol{\nu}}\end{array} +$$ + +These coefficients and the generalized force are only functions of time $t$ and the displacements ${\bf q}(t)$ . The first term of (1.7) describes the fundamental inertia force due to accelerations $\ddot{\mathbf{q}}$ with a symmetric mass matrix $m_{i j}=m_{j i}$ . The second term describes the gyroscopic forces from substructures where the position vector has a explicit dependency of time $\partial\mathbf{r}/\partial t\neq0$ , for example if the kinematic formulation is based on a drivetrain that rotates at a given average speed and the drivetrain dynamics is described by variations around this speed. The third term describes the centrifugal forces on substructures rotating with a prescribed speed, or the forces due to other explicitly defined accelerations such as earthquakes or the moving base when this structural model is modular coupled to another dynamic model of a foundation or floater. The fourth term can describe the similar Coriolis and acceleration forces as the second and third terms. For example let the generalized coordinate $q_{k}$ be the absolute rotation angle of the drivetrain then this kinematic formulation can be changed to one with a prescribed averaged rotational speed by the substitution $q_{k}\,=\,\Omega t+\delta q_{k}$ , where $\Omega$ is the prescribed average speed and $\delta q_{k}$ is the new generalized coordinate. The part term $h_{i j k}\dot{q}_{k}=h_{i j k}\Omega$ will have a component equal to the coefficient $g_{i j}$ in case of a prescribed constant speed bearing for $j\neq k$ , and the entire term $h_{i j k}\dot{q}_{j}\dot{q}_{k}=h_{i j k}\Omega^{2}$ for $j=k$ will have a component equal to the acceleration force $F_{c,i}$ . + +In case of a prescribed rotation of the rotor, the acceleration forces $F_{c,i}$ given by (1.8d) are centrifugal forces that stiffen the blades. To include this centrifugal stiffness in the calculation of frequencies or in the iteration steps of a time integration, we can compute the centrifugal stiffness matrix as the Jacobian of this vector function as + +$$ +k_{c,i j}=\int_{\mathcal{V}}\rho\left(\frac{\partial^{2}\pmb{r}^{T}}{\partial q_{i}\partial q_{j}}\frac{\partial^{2}\pmb{r}}{\partial t^{2}}+\frac{\partial\pmb{r}^{T}}{\partial q_{i}}\frac{\partial^{3}\pmb{r}}{\partial t^{2}\partial q_{j}}\right)d\mathcal{V} +$$ + +which is not a symmetric matrix due to the second term. + +The fifth term of (1.7) describes the conservative forces such as gravity and elastic forces which can be defined by the potential energy $V$ . The force function given by the derivative of the scalar potential energy function depends on the kinematic formulation and the applied theory for elasticity. The sixth term of (1.7) describes the purely dissipative forces from structural damping mechanisms (e.g. material and friction), which we will model using a modal damping method described in Section... . The last term of (1.7) describes generalized non­conservative forces acting on the structure which are discussed in Chapter... . + +# 1.1.1 Generic topology of structure + +The turbine structure is divided into a number of articulated substructures where one substructure may be connected to one or more substructures through its connection points. Let $b$ be the index of a substructure then we herein assume that it is connected to substructure $b-1$ which is connected to substructure $b-2$ and so on until substructure 0. In practice, this numbering will not be continuous (e.g. the blades will each have their own number but connected to the same hub which can have one number), but we can always create an intermediate index list for each substructure in which the numbering is continuous. + +The position vector $\boldsymbol{\mathsf{r}}_{b}$ of a particle on substructure number $b$ is written as + +$$ +\begin{array}{r l}&{{\bf{r}}_{b}\left(t,{\bf{q}};x,y,z\right)=\!{\bf{S}}_{0}{\bf{B}}_{0}\left({\bf{r}}_{0}+{\bf{R}}_{0}{\bf{S}}_{0}^{T}{\bf{S}}_{1}{\bf{B}}_{1}\left({\bf{r}}_{1}+{\bf{R}}_{1}{\bf{S}}_{1}^{T}{\bf{S}}_{2}{\bf{B}}_{2}\left({\bf{r}}_{2}\right.\right.}\\ &{\quad\quad\quad\quad\quad\quad\left.+{\bf{R}}_{2}{\bf{S}}_{2}^{T}{\bf{S}}_{3}{\bf{B}}_{3}\left({\bf{r}}_{3}+\cdot\cdot+{\bf{R}}_{b-1}{\bf{S}}_{b-1}^{T}{\bf{S}}_{b}{\bf{B}}_{b}\left.{\bf{r}}_{1,b}\left({\bf{q}}_{b};x,y,z\right)\right.\right.}\end{array} +$$ + +where $t$ is time, $x,y$ and $z$ are the coordinates of the particle in the moving substructure frame and $\mathbf{q}$ is the DOF vector. Note that $\mathbf{q}_{b}$ is a vector containing the subset of $\mathbf{q}$ with the DOFs of the substructure $b$ . The matrices $\pmb{\mathrm{s}}_{b}$ are the constant rotation matrix for the coordinate systems of the substructure $b$ described in the ground fixed inertia frame, $\mathbf{B}_{b}$ are the rotation matrix for a possible bearing in the base1 of substructure $b$ , and $\mathbf{R}_{b}$ are the rotation matrix at the connection point on substructure $b$ to which the next substructure $b+1$ is connected. The vectors $\boldsymbol{\mathsf{r}}_{b}$ describe the positions of these connection points in the coordinates of substructure $b$ , and the vector $\boldsymbol{\mathsf{r}}_{1,b}$ is internal deformation vector function describing the position of all particles in the entire substructure $b$ in its local frame of reference. This internal deformation is given by the DOFs $\mathbf{q}_{b}$ . The vector $\boldsymbol{\mathsf{r}}_{b}$ and the matrices $\mathbf{B}_{b}$ and $\mathbf{R}_{b}$ are functions of the DOFs $\mathbf{q}_{b}$ . For prescribed bearings, e.g. constant rotation speed, the matrix $\mathbf{B}_{b}$ may also an explicit function of time $t$ . Without prescribed bearings, there are no explicit functions of time in the position vectors for any substructure, whereby the gyroscopic matrix and centrifugal force vector vanish. For a free bearing with a DOF describing the rotation of the drivetrain, the gyroscopic and centrifugal inertia forces are included in the nonlinear Coriolis term of (1.7) given by the coefficients in (1.8c). + +The position vector (1.10) can be rewritten in condensed form as + +$$ +\mathbf{r}_{b}\left(t,\mathbf{q};x,y,z\right)=\sum_{k=0}^{b-1}\left(\prod_{l=0}^{k-1}\mathbf{O}_{l}\right)\mathbf{S}_{k}\mathbf{B}_{k}\mathbf{r}_{k}+\left(\prod_{l=0}^{b-1}\mathbf{O}_{l}\right)\mathbf{S}_{b}\mathbf{B}_{b}\mathbf{r}_{1,b}\left(\mathbf{q}_{b};x,y,z\right) +$$ + +where $\begin{array}{r}{\pmb{0}_{l}=\pmb{\mathsf{S}}_{l}\pmb{\mathsf{B}}_{l}\pmb{\mathsf{R}}_{l}\pmb{\mathsf{S}}_{l}^{T}}\end{array}$ is the rotation matrix of substructure $l$ described in the ground­fixed frame, including both the deformations of the substructure and the rotations of the bearing at its base. For the following derivations, the position vector (1.11) is written in the further condensed form as + +$$ +\begin{array}{r}{\mathbf{r}_{b}\left(t,\mathbf{q};x,y,z\right)=\mathbf{r}_{0,b}\left(t,\mathbf{q}\right)+\mathbf{R}_{0,b}\left(t,\mathbf{q}\right)\mathbf{r}_{1,b}\left(\mathbf{q}_{b};x,y,z\right)}\end{array} +$$ + +where the translations and rotations of the substructure base are described by + +$$ +\mathbf{r}_{0,b}=\sum_{k=0}^{b-1}\mathbf{R}_{0,k}\mathbf{r}_{k}=\mathbf{r}_{0,b-1}+\mathbf{R}_{0,b-1}\mathbf{r}_{b-1}\quad{\mathrm{and}}\quad\mathbf{R}_{0,b}=\left(\prod_{l=0}^{b-1}\mathbf{O}_{l}\right)\mathbf{S}_{b}\mathbf{B}_{b}=\mathbf{R}_{0,b-1}\mathbf{R}_{b-1}\mathbf{S}_{b-1}^{T}\mathbf{S}_{b}\mathbf{B}_{b} +$$ + +where $\mathsf{r}_{k}$ is the local position vector in the frame of substructure $k$ to which substructure $k+1$ is connected; this vector is independent of time because there are no bearings involved. We use the index “1” for the local substructure deformation given by $\boldsymbol{\mathsf{r}}_{1,b}$ and $^{\ast}0^{\ast}$ for the deformation and orientation of its base given by $\boldsymbol{\mathsf{r}}_{0,b}$ and $\pmb{\mathsf{R}}_{0,b}$ . Note that these position vector and orientation matrix are functions of DOFs of substructures supporting the substructure $b$ , but not of DOFs of other substructures or substructure $b$ itself. + +In the subsequent derivations, we will need the time derivatives: + +$$ +\begin{array}{r l}&{\tilde{F}_{\mathrm{0,A}}=\displaystyle\sum_{k=1}^{k-1}\mathbb{E}_{\boldsymbol{u},\boldsymbol{u}^{k}}}\\ &{\tilde{\mathbf{R}}_{\mathrm{0,A}}=\displaystyle\sum_{k=1}^{k-1}\left(\prod_{i=0}^{n-1}\mathbf{\sigma}_{0}\right)\le\mathbf{1}_{\mathcal{B}_{i}\setminus\mathbb{R}_{i}}\mathbf{1}_{\mathcal{S}_{i}^{\top}}\left(\prod_{i=1}^{k-1}\mathbf{\sigma}_{0}\right)\le\mathbf{1}_{\mathcal{B}_{i}\setminus+}\left(\prod_{i=0}^{n-1}\mathbf{\sigma}_{0}\right)\le\mathbf{1}_{\mathcal{B}_{i}},}\\ &{\tilde{F}_{\mathrm{0,A}}=\displaystyle\sum_{k=0}^{k-1}\mathbb{E}_{\boldsymbol{u},\boldsymbol{u}^{k}}}\\ &{\tilde{F}_{\mathrm{0,A}}=\displaystyle\sum_{k=1}^{k-1}\mathbb{E}_{\boldsymbol{u},\boldsymbol{u}^{k}}}\\ &{\tilde{\mathbf{R}}_{\mathrm{0,B}}=\displaystyle\sum_{k=0}^{k-1}\left(\sum_{i=0}^{k-1}\bigg(\prod_{i=1}^{n}\mathbf{\sigma}_{i}\bigg)\le\mathbf{1}_{\mathcal{B}_{i}\setminus\mathbb{R}_{i}}\mathbf{1}_{\mathcal{S}_{i}^{\top}}\left(\prod_{i=1}^{k-1}\mathbf{\sigma}_{i}\right)\right)\le\mathbf{1}_{\mathcal{B}_{i}\setminus\mathbb{R}_{i}}\mathbf{1}_{\mathcal{S}_{i}^{\top}}\left(\prod_{i=1}^{k-1}\mathbf{\sigma}_{i}\right)\le\mathbf{1}_{\mathcal{B}_{i}},}\\ &{\qquad+\displaystyle\sum_{k=1}^{k-1}\left(\prod_{i=0}^{k-1}\mathbf{\sigma}_{i}\right)\le\mathbf{1}_{\mathcal{B}_{i}\setminus\mathbb{R}_{i}}\mathbf{1}_{\mathcal{S}_{i}^{\top}}\left(\begin{array}{c c c}{\displaystyle\sum_{i=1}^{k-1}\left(\prod_{i=1}^{n}\mathbf{\sigma}_{i}\right)}&{\mathbf{S}_{i}\mathbf{1}_{\mathcal{B}_{i}}\mathbf{1}_{\mathcal{S}_{i}^{\top}} +$$ + +where we have used that the bearing matrices $\mathbf{B}_{k}$ represent the only explicit time dependency. Note that some bearings have a constant angle or are function of a free bearing DOF, whereby $\dot{\ B_{k}}=\ddot{\ B}_{k}=\pmb{0}$ . In most applications, there will only be a single bearing with an explicit time dependency in terms of a constant speed of $\Omega$ . Let this single bearing be for substructure number $k$ then the time derivatives (1.14) reduce to + +$$ +\begin{array}{r l}&{\dot{\mathbf{r}}_{0,b}=\displaystyle\sum_{l=k}^{b-1}\dot{\mathbf{R}}_{0,l}\mathbf{r}_{l}}\\ &{\dot{\mathbf{R}}_{0,b}=\Omega\left(\prod_{l=0}^{k-1}\mathbf{O}_{l}\right)\,\mathbf{S}_{k}\mathbf{B}_{k}\mathbf{N}_{k}\mathbf{R}_{k}\mathbf{S}_{k}^{T}\,\left(\displaystyle\prod_{l=k+1}^{b-1}\mathbf{O}_{l}\right)\mathbf{S}_{b}\mathbf{B}_{b}=\Omega\mathbf{R}_{0,k}\mathbf{N}_{k}\mathbf{R}_{0,k}^{T}\mathbf{R}_{0,b}}\\ &{\ddot{\mathbf{r}}_{0,b}=\displaystyle\sum_{l=k}^{b-1}\dot{\mathbf{R}}_{0,l}\mathbf{r}_{l}}\\ &{\ddot{\mathbf{R}}_{0,b}=\Omega^{2}\left(\prod_{l=0}^{k-1}\mathbf{O}_{l}\right)\,\mathbf{S}_{k}\mathbf{B}_{k}\mathbf{N}_{k}^{2}\mathbf{R}_{k}\mathbf{S}_{k}^{T}\,\left(\displaystyle\prod_{l=k+1}^{b-1}\mathbf{O}_{l}\right)\,\mathbf{S}_{b}\mathbf{B}_{b}=\Omega^{2}\mathbf{R}_{0,k}\mathbf{N}_{k}^{2}\mathbf{R}_{0,k}^{T}\mathbf{R}_{0,l}}\end{array} +$$ + +where $\mathbf{N}_{k}$ is the skew symmetric matrix (E.3) defined by the unit­vectors of the bearing axis. In this case, the following products that occur in the subsequent derivations can be written as + +$$ +\begin{array}{r l}&{\mathbf{r}_{0,b}^{T}\mathbf{\Phi}_{0,b}^{\dagger}=}\\ &{\mathbf{R}_{0,b}^{T}\mathbf{\Phi}_{0,b}^{\dagger}=\Omega\mathbf{R}_{0,b}^{T}\mathbf{R}_{0,i}\mathbf{N}_{k}^{\mathbf{R}}\mathbf{R}_{0,b}^{T}\sum_{l=k}^{b-1}\mathbf{R}_{0,l}\mathbf{r}_{l}=\Omega\left(\mathbf{R}_{0,k}^{T}\mathbf{R}_{0,b}\right)^{T}\mathbf{N}_{k}\mathbf{R}_{0,k}^{T}\left(\mathbf{r}_{0,b}-\mathbf{r}_{0,k}\right)}\\ &{\mathbf{R}_{0,b}^{T}\mathbf{\Phi}_{0,b}^{\dagger}=\Omega\left(\mathbf{R}_{0,k}^{T}\mathbf{R}_{0,b}\right)^{T}\mathbf{N}_{k}\mathbf{\Phi}_{0,b}^{T}\mathbf{R}_{0,b}}\\ &{\mathbf{r}_{0,b}^{T}\mathbf{\Phi}_{0,b}^{\dagger}=}\\ &{\mathbf{R}_{0,b}^{T}\mathbf{\Phi}_{0,b}^{\dagger}=\Omega^{2}\mathbf{R}_{0,b}^{T}\mathbf{R}_{0,b}^{2}\sum_{l=k}^{b-1}\mathbf{R}_{0,l}\mathbf{r}_{l}=\Omega^{2}\left(\mathbf{R}_{0,k}^{T}\mathbf{R}_{0,b}\right)^{T}\mathbf{N}_{k}^{2}\mathbf{R}_{0,k}^{T}\left(\mathbf{r}_{0,b}-\mathbf{r}_{0,k}\right)}\\ &{\mathbf{R}_{0,b}^{T}\mathbf{\Phi}_{0,b}^{\dagger}=\Omega^{2}\left(\mathbf{R}_{0,b}^{T}\mathbf{R}_{0,b}\right)^{T}\mathbf{R}_{k}^{2}\mathbf{R}_{0,k}^{T}\mathbf{R}_{0,b}}\end{array} +$$ + +which are time independent scalers, vectors and matrices because + +$$ +\mathbf{\mathfrak{I}}_{0,k}^{T}\mathbf{R}_{0,b}=\mathbf{B}_{k}^{T}\mathbf{S}_{k}^{T}\left(\prod_{l=k-1}^{0}\mathbf{O}_{l}^{T}\right)\left(\prod_{l=0}^{b-1}\mathbf{O}_{l}\right)\mathbf{S}_{b}\mathbf{B}_{b}=\mathbf{B}_{k}^{T}\mathbf{S}_{k}^{T}\left(\prod_{l=k}^{b-1}\mathbf{O}_{l}\right)\mathbf{S}_{b}\mathbf{B}_{b}=\mathbf{R}_{k}\mathbf{S}_{k}^{T}\left(\prod_{l=k+1}^{b-1}\mathbf{O}_{l}\right)\mathbf{S}_{b}\mathbf{B}_{b} +$$ + +if $b\geq k$ + +We also need the DOF derivatives: + +$$ +\begin{array}{r l}{\displaystyle\mathbf{r}_{0,b,q_{i}}=\sum_{k=0}^{b-1}\big(\mathbf{R}_{0,k,q_{i}}\mathbf{r}_{k}+\mathbf{R}_{0,k}\mathbf{r}_{k,q_{i}}\big)}&{}\\ {\displaystyle\mathbf{R}_{0,b,q_{i}}=}&{}\\ {\displaystyle\mathbf{r}_{0,b,q_{i},q_{j}}=}&{}\\ {\displaystyle\mathbf{R}_{0,b,q_{i},q_{j}}=}&{}\end{array} +$$ + +# 1.1.2 Inertia forces on and from a substructure + +The inertia forces from a substructure is “felt” by all substructures supporting it, i.e., there may only be entries in rows and columns of the mass, gyroscopic, and centrifugal stiffness matrices for DOFs of the supporting substructures that appear in the vector and matrix functions $\boldsymbol{\mathsf{r}}_{0,b}$ and $\mathbf{R}_{0,b}$ . We subdivide the mass, gyroscopic, and centrifugal stiffness matrices into $_{2\times2}$ block matrices, e.g. the mass matrix as + +$$ +\boldsymbol{\mathsf{M}}=\left[\begin{array}{l l}{\boldsymbol{\mathsf{M}}^{00}}&{\boldsymbol{\mathsf{M}}^{01}}\\ {\boldsymbol{\mathsf{M}}^{10}}&{\boldsymbol{\mathsf{M}}^{11}}\end{array}\right] +$$ + +where the diagonal matrices $\mathbb{M}^{00}$ and $\mathbb{M}^{11}$ are the mass matrix contributions for the DOFs of the supporting sub­ structures and the DOFs of the substructure, respectively, and $\mathbb{M}^{10}$ and $\pmb{\mathbb{M}}^{01}$ are the mass matrix contributions that coupled these DOFs. Remember that $\mathbb{M}^{10}=\mathbb{M}^{01^{T}}$ because the mass matrix is symmetric; the gyroscopic and centrifugal stiffness matrices are not symmetric. + +Similar, we subdivide the nonlinear three­dimensional gyroscopic matrix into two $_{2\times2}$ block matrices: + +$$ +\mathsf{H}_{i}^{0}=\left[\begin{array}{l l}{\mathsf{H}^{000}}&{\mathsf{H}^{001}}\\ {\mathsf{H}^{010}}&{\mathsf{H}^{011}}\end{array}\right]\ \mathsf{a n d}\ \mathsf{H}_{i}^{1}=\left[\begin{array}{l l}{\mathsf{H}^{100}}&{\mathsf{H}^{101}}\\ {\mathsf{H}^{110}}&{\mathsf{H}^{111}}\end{array}\right] +$$ + +where $\mathsf{H}_{i}^{0}$ contains the nonlinear gyroscopic matrix elements for DOF $q_{i}$ on a supporting substructure, and $\mathsf{H}_{i}^{1}$ contains the elements for DOF $q_{i}$ on the substructure itself. Thus, the element $h_{i j k}^{000}$ describe the coefficient of the gyroscopic force from velocities ${\dot{q}}_{j}$ and $\dot{q}_{k}$ of one or two different supporting substructures on the DOF $q_{i}$ of the same or different supporting substructure. + +In the following, the substructure­index $b$ are omitted for brevity. + +# Mass matrix + +Inserting (1.12) into (1.8a), expanding and sorting the terms into the $_{2\times2}$ block matrix form yield the mass matrix elements + +$$ +\begin{array}{r}{m_{i j}^{00}=\int_{\mathcal{V}}\rho\left(\mathbf{r}_{0,q_{i}}^{T}\mathbf{r}_{0,q_{j}}+\left(\mathbf{r}_{0,q_{i}}^{T}\mathbf{R}_{0,q_{j}}+\mathbf{r}_{0,q_{j}}^{T}\mathbf{R}_{0,q_{i}}\right)\mathbf{r}_{1}+\mathbf{r}_{1}^{T}\mathbf{R}_{0,q_{i}}^{T}\mathbf{R}_{0,q_{j}}\mathbf{r}_{1}\right)d\mathcal{V}}\\ {m_{i j}^{01}=m_{j i}^{10}=\int_{\mathcal{V}}\rho\left(\mathbf{r}_{0,q_{i}}^{T}\mathbf{R}_{0}\mathbf{r}_{1,q_{j}}+\mathbf{r}_{0,q_{j}}^{T}\mathbf{R}_{0}\mathbf{r}_{1,q_{i}}+\mathbf{r}_{1}^{T}\mathbf{R}_{0,q_{i}}^{T}\mathbf{R}_{0}\mathbf{r}_{1,q_{j}}+\mathbf{r}_{1}^{T}\mathbf{R}_{0,q_{j}}^{T}\mathbf{R}_{0}\mathbf{r}_{1,q_{i}}\right)d\mathcal{V}}\\ {m_{i j}^{11}=\int_{\mathcal{V}}\rho\left(\mathbf{r}_{1,q_{i}}^{T}\mathbf{r}_{1,q_{j}}\right)d\mathcal{V}}\end{array} +$$ + +where we have used the notation $()_{,q_{i}}\equiv\partial/\partial q_{i}$ for the first derivatives with respect to DOF $q_{i}$ . It is possible to isolate the integration over the entire substructure volume to the local deformation vector and its derivatives by using the formula (E.7) for the matrix dot product $\equiv\mathsf{A}:\mathsf{B}$ . The mass matrix elements can thereby be written as + +$$ +\begin{array}{r l}&{\quad m_{i j}^{00}=\!\!\Gamma_{0,q_{i}}^{T}\mathbf{r}_{0,q_{j}}\!\int_{\gamma}\rho d\gamma+\left(\mathbf{r}_{0,q_{i}}^{T}\mathbf{R}_{0,q_{j}}+\mathbf{r}_{0,q_{j}}^{T}\mathbf{R}_{0,q_{i}}\right)\int_{\gamma}\rho\mathbf{r}_{1}d\gamma+\left(\mathbf{R}_{0,q_{i}}^{T}\mathbf{R}_{0,q_{j}}\right):\int_{\gamma}\rho\mathbf{r}_{1}\mathbf{r}_{1}^{T}d\gamma}\\ &{\quad m_{i j}^{01}=m_{j i}^{10}=\!\!\Gamma_{0,q_{i}}^{T}\mathbf{R}_{0}\!\int_{\gamma}\rho\mathbf{r}_{1,q_{j}}d\gamma+\mathbf{r}_{0,q_{j}}^{T}\mathbf{R}_{0}\!\int_{\gamma}\rho\mathbf{r}_{1,q_{i}}d\gamma}\\ &{\quad\quad\quad\quad\quad+\left(\mathbf{R}_{0,q_{i}}^{T}\mathbf{R}_{0}\right):\int_{\gamma}\rho\mathbf{r}_{1,q_{j}}\mathbf{r}_{1}^{T}d\gamma+\left(\mathbf{R}_{0,q_{j}}^{T}\mathbf{R}_{0}\right):\int_{\gamma}\rho\mathbf{r}_{1,q_{i}}\mathbf{r}_{1}^{T}d\gamma}\\ &{\quad\quad\quad\quad\quad\quad\quad\quad m_{i j}^{11}=\!\!\int_{\gamma}\rho\left(\mathbf{r}_{1,q_{i}}^{T}\mathbf{r}_{1,q_{j}}\right)d\gamma}\end{array} +$$ + +where the blue colored integrals of these expressions can computed in the code object of the substructure, independent of its base motion. Some of these integrals over the volume of the substructure have physical meanings as + +$$ +\begin{array}{r l}&{\quad\quad\displaystyle\int_{\mathcal{V}}\rho d\mathcal{V}=M}\\ &{\quad\quad\displaystyle\int_{\mathcal{V}}\rho\,\mathbf{r}_{1}d\mathcal{V}=M\left\{\begin{array}{l}{x_{c g}}\\ {y_{c g}}\\ {z_{c g}}\end{array}\right\}=M\mathbf{r}_{c g}}\\ &{\quad\displaystyle\int_{\mathcal{V}}\rho\,\mathbf{r}_{1}\mathbf{r}_{1}^{T}d\mathcal{V}=\left[\begin{array}{l l l}{i_{x x}}&{i_{x y}}&{i_{x z}}\\ {i_{x y}}&{i_{y y}}&{i_{y z}}\\ {i_{x z}}&{i_{y z}}&{i_{z z}}\end{array}\right]=\mathbf{l}_{b a s e}}\end{array} +$$ + +where $M$ is the total mass of the substructure, $\mathbf{r}_{c g}=\{x_{c g},y_{c g},z_{c g}\}^{T}$ is the center of gravity of the substructure measured from its base in the ground­fixed inertia frame, and $\mid_{b a s e}$ is a matrix related to the rotational inertia of + +the substructure about its base given by the integrals defined as + +$$ +i_{\alpha\beta}=\int_{\mathcal{V}}\rho\,r_{1,\alpha}r_{1,\beta}d\mathcal{V} +$$ + +where the subscripts $(\alpha,\beta)\in\left\{x,y,z\right\}$ denote the coordinate component of the position vector. These integrals can be used to compute the moments and products of inertia [1] for the substructure about its base defined in the ground­fixed frame as + +$$ +\begin{array}{c}{I_{x x}=i_{y y}+i_{z z}\ ,\ I_{y y}=i_{x x}+i_{z z}\ ,\ I_{z z}=i_{x x}+i_{y y}\ ,}\\ {I_{x y}=I_{y x}=-i_{x y}\ ,\ I_{x z}=I_{z x}=-i_{x z}\ ,\mathsf{a n d}\ I_{y z}=I_{z y}=-i_{y z}}\end{array} +$$ + +The remaining integrals over the substructure volume of (1.22) involve the derivatives of the local position vector with respect to the DOFs of the substructure. One is the derivative of the center of gravity position: + +$$ +\int_{\mathcal{V}}\rho\,\mathbf{r}_{1,q_{i}}d\mathcal{V}=M\mathbf{r}_{c g,q_{i}} +$$ + +Another volume integral is the “asymmetric half” of the derivative of the symmetric inertia matrix + +$$ +\mathbf{l}_{b a s e,q_{i}}=\int_{\mathcal{V}}\rho\,\mathbf{r}_{1,q_{i}}\mathbf{r}_{1}^{T}d\mathcal{V}+\int_{\mathcal{V}}\rho\,\mathbf{r}_{1}\mathbf{r}_{1,q_{i}}^{T}d\mathcal{V}=\int_{\mathcal{V}}\rho\,\mathbf{r}_{1,q_{i}}\mathbf{r}_{1}^{T}d\mathcal{V}+\left(\int_{\mathcal{V}}\rho\,\mathbf{r}_{1,q_{i}}\mathbf{r}_{1}^{T}d\mathcal{V}\right)^{T}=\mathbf{A}_{b a s e,i}+\mathbf{A}_{b a s e,i}^{T} +$$ + +where the notation $\begin{array}{r}{\pmb{\mathsf{A}}_{b a s e,i}\;\equiv\;\int_{\mathcal{V}}\!\rho\,\pmb{\mathsf{r}}_{1,q_{i}}\pmb{\mathsf{r}}_{1}^{T}d\mathcal{V}}\end{array}$ is introduced for this volume integral. Finally, the last integral in (1.22c) defined the entries of the local mass matrix of the substructure $\pmb{\mathbb{M}}^{11}$ independent of the base motion. + +# Gyroscopic matrix + +Inserting (1.12) into (1.8b), expanding and sorting the terms into the $_{2\times2}$ block matrix form yield the gyroscopic matrix elements + +$$ +\begin{array}{l}{{g_{i j}^{00}=2\displaystyle\int_{\mathcal{V}}\rho\left(\mathbf{r}_{0,q_{i}}^{T}\dot{\mathbf{r}}_{0,q_{j}}+\mathbf{r}_{0,q_{i}}^{T}\dot{\mathbf{R}}_{0,q_{j}}\mathbf{r}_{1}+\dot{\mathbf{r}}_{0,q_{j}}^{T}\mathbf{R}_{0,q_{i}}\mathbf{r}_{1}+\mathbf{r}_{1}^{T}\mathbf{R}_{0,q_{i}}^{T}\dot{\mathbf{R}}_{0,q_{j}}\mathbf{r}_{1}\right)d\mathcal{V}}}\\ {{g_{i j}^{01}=2\displaystyle\int_{\mathcal{V}}\rho\left(\dot{\mathbf{r}}_{0,q_{j}}^{T}\mathbf{R}_{0}\mathbf{r}_{1,q_{i}}+\mathbf{r}_{0,q_{i}}^{T}\dot{\mathbf{R}}_{0}\mathbf{r}_{1,q_{j}}+\mathbf{r}_{1}^{T}\mathbf{R}_{0,q_{i}}^{T}\dot{\mathbf{R}}_{0}\mathbf{r}_{1,q_{j}}+\mathbf{r}_{1}^{T}\dot{\mathbf{R}}_{0,q_{j}}^{T}\mathbf{R}_{0}\mathbf{r}_{1,q_{i}}\right)d\mathcal{V}}}\\ {{g_{i j}^{11}=2\displaystyle\int_{\mathcal{V}}\rho\mathbf{r}_{1,q_{i}}^{T}\mathbf{R}_{0}^{T}\dot{\mathbf{R}}_{0}\mathbf{r}_{1,q_{j}}d\mathcal{V}}}\end{array} +$$ + +and $g_{i j}^{10}$ is simply given by switching the indices $i$ and $j$ in the expression (1.28b) for $g_{i j}^{01}$ . We isolate the integration over the entire substructure volume to the local deformation vector and its derivatives by using (E.7), whereby the gyroscopic matrix elements are written as + +$$ +\begin{array}{r l}&{g_{i j}^{00}=2M\mathbf{r}_{0,q_{i}}^{T}\dot{\mathbf{t}}_{0,q_{j}}+2\left(\mathbf{r}_{0,q_{i}}^{T}\dot{\mathbf{R}}_{0,q_{j}}+\dot{\mathbf{r}}_{0,q_{j}}^{T}\mathbf{R}_{0,q_{i}}\right)M\mathbf{r}_{c g}+2\left(\mathbf{R}_{0,q_{i}}^{T}\dot{\mathbf{R}}_{0,q_{j}}\right):\mathbf{l}_{b a s e}}\\ &{g_{i j}^{01}=2\dot{\mathbf{r}}_{0,q_{j}}^{T}\mathbf{R}_{0}M\mathbf{r}_{c g,q_{i}}+2\mathbf{r}_{0,q_{i}}^{T}\dot{\mathbf{R}}_{0}M\mathbf{r}_{c g,q_{j}}+2\left(\mathbf{R}_{0,q_{i}}^{T}\dot{\mathbf{R}}_{0}\right):\mathbf{A}_{b a s e,j}+2\left(\dot{\mathbf{R}}_{0,q_{j}}^{T}\mathbf{R}_{0}\right):\mathbf{A}_{b a s e,i}}\\ &{g_{i j}^{11}=2\left(\mathbf{R}_{0}^{T}\dot{\mathbf{R}}_{0}\right):\mathsf{A}_{b a s e,1,i j}}\end{array} +$$ + +where $\begin{array}{r}{\mathsf{A}_{b a s e,1,i j}\;\equiv\;\int_{\mathcal{V}}\!\rho\,\mathsf{r}_{1,q_{j}}\mathsf{r}_{1,q_{i}}^{T}d\mathcal{V}}\end{array}$ is introduced for the volume integral over the substructure of the matrix defined by the first derivatives of the local deformation vector with respect to the local DOFs. + +# Nonlinear gyroscopic matrix + +Inserting (1.12) into (1.8c), expanding and sorting the terms into the two $_{2\times2}$ block matrix forms (1.20) yield the nonlinear gyroscopic matrix elements + +$$ +\begin{array}{r l}&{\tilde{h}_{i j k}^{\mathrm{OO}}=\int_{\gamma}\left(\tau_{\mathrm{O}_{i}}^{T}\mathbf{R}_{\mathrm{O}_{j},\mathrm{f}_{i}+1}+\tau_{\mathrm{O}_{j}}^{T}\mathbf{R}_{\mathrm{O}_{i},\mathrm{f}_{i}+1}+\tau_{\mathrm{O}_{j},\mathrm{f}_{i}+1}^{T}\mathbf{R}_{\mathrm{O}_{j},\mathrm{f}_{i}+1}+\tau_{\mathrm{P}_{i}}^{T}\mathbf{R}_{\mathrm{O}_{j},\mathrm{f}_{i}+1}^{T}\right)d\gamma}\\ &{\tilde{h}_{i j k}^{\mathrm{O}}=\int_{\gamma}\rho\left(\tau_{\mathrm{O}_{i}}^{T}\mathbf{R}_{\mathrm{O}_{j},\mathrm{f}_{i}+1}+\tau_{\mathrm{P}_{i}}^{T}\mathbf{R}_{\mathrm{O}_{j},\mathrm{f}_{i}+1}^{T}\right)d\gamma}\\ &{\tilde{h}_{i j k}^{\mathrm{O}}=\int_{\gamma}\rho\left(\tau_{\mathrm{O}_{i}}^{T}\mathbf{R}_{\mathrm{O}_{j},\mathrm{f}_{i}+1}+\tau_{\mathrm{P}_{i}}^{T}\mathbf{R}_{\mathrm{O}_{j},\mathrm{f}_{i}+1}^{T}\right)d\gamma}\\ &{\tilde{h}_{i j k}^{\mathrm{O}}=\int_{\gamma}\rho\left(\tau_{\mathrm{O}_{i}}^{T}\mathbf{R}_{\mathrm{O}_{j},\mathrm{f}_{i}+1}+\tau_{\mathrm{P}_{i}}^{T}\mathbf{R}_{\mathrm{O}_{j},\mathrm{f}_{i}+1}^{T}\right)d\gamma}\\ &{\tilde{h}_{i j k}^{\mathrm{O}}=\int_{\gamma}\rho\left(\tau_{\mathrm{O}_{i},\mathrm{f}_{i}+1}^{T}\mathbf{R}_{\mathrm{O}_{j},\mathrm{f}_{i}+1}+\tau_{\mathrm{P}_{i}}^{T}\mathbf{R}_{\mathrm{O}_{j},\mathrm{f}_{i}+1}^{T}\right)d\gamma}\\ &{\tilde{h}_{i j k}^{\mathrm{IO}}=\int_{\gamma}\rho\left(\tau_{\mathrm{O}_{i},\mathrm{f}_{i}+1}^{T}\mathbf{R}_{ +$$ + +where we have used the notation $()_{,q_{i},q_{j}}\equiv\partial^{2}/\partial q_{i}\partial q_{j}$ for the second derivatives with respect to DOFs $q_{i}$ and $q_{j}$ . Again, the integration over the entire substructure volume is isolated to the local deformation vector and its derivatives by using (E.7), whereby the nonlinear gyroscopic matrix elements are written as + +$$ +\begin{array}{r l}&{h_{i j k}^{00}=M_{0,q_{i}}^{T}\mathbf{\Delta}\mathbf{f}_{0,q_{j},q_{k}}+\left(\mathbf{r}_{0,q_{i}}^{T}\mathbf{R}_{0,q_{j},q_{k}}+\mathbf{r}_{0,q_{i},q_{k}}^{T}\mathbf{R}_{0,q_{i}}\right)M_{\mathbf{r},q}+\left(\mathbf{R}_{0,q_{i}}^{T}\mathbf{R}_{0,q_{j},q_{k}}\right):\mathbf{\Delta}_{\mathbf{b}_{\mathrm{asc}}}}\\ &{h_{i j k}^{00}=\mathbf{r}_{0,q_{i}}^{T}\mathbf{R}_{0,q_{j}}M_{\mathbf{r},q_{k}}+\left(\mathbf{R}_{0,q_{i}}^{T}\mathbf{R}_{0,q_{j}}\right):\mathbf{\Delta}\mathbf{\hat{A}}_{\mathbf{b}_{\mathrm{asc}},k}}\\ &{h_{i j k}^{01}=\mathbf{r}_{0,q_{i}}^{T}\mathbf{R}_{0,q_{k}}M_{\mathbf{r},q_{j},q_{k}}+\left(\mathbf{R}_{0,q_{i}}^{T}\mathbf{R}_{0,q_{k}}\right):\mathbf{\Delta}\mathbf{\hat{A}}_{\mathbf{b}_{\mathrm{asc}},j}}\\ &{h_{i j k}^{01}=\mathbf{r}_{0,q_{i}}^{T}\mathbf{R}_{0}M_{\mathbf{r},q_{j},q_{k}}+\left(\mathbf{R}_{0,q_{i}}^{T}\mathbf{R}_{0}\right):\mathbf{\Delta}\mathbf{\hat{A}}_{\mathbf{basc},j,j}}\\ &{h_{i j k}^{100}=\mathbf{r}_{0,q_{i},q_{k}}^{T}\mathbf{R}_{0}M\mathbf{r}_{\mathbf{c},q_{i}}+\left(\mathbf{R}_{0,q_{i},q_{k}}^{T}\mathbf{R}_{0}\right):\mathbf{\Delta}\mathbf{\hat{A}}_{\mathbf{basc},i,j}}\\ &{h_{i j k}^{101}=\left(\mathbf{R}_{0}^{T}\mathbf{R}_{0,q_{j}}\right):\mathbf{\Delta}\mathbf{\hat{A}}_{\mathbf{basc},1,i}}\\ &{h_{i j k}^{11}=\left( +$$ + +where $\begin{array}{r}{\mathsf{A}_{b a s e,2,i j}\;\equiv\;\int_{\mathcal{V}}\rho\,\mathsf{r}_{1,q_{i},q_{j}}\,\mathsf{r}_{1}^{T}d\mathcal{V}}\end{array}$ is introduced for the volume integral over the substructure of the matrix defined by the local deformation vector and its second derivatives with respect to the local DOFs. + +# Centrifugal force vector + +Inserting (1.12) into (1.8d) and expanding lead to these acceleration (centrifugal) forces on the supporting DOFs and the local substructure DOFs: + +$$ +\begin{array}{l}{{\displaystyle F_{c,i}^{0}=\int_{\mathcal{V}}\rho\left({\bf r}_{0,q_{i}}^{T}\ddot{\bf r}_{0}+\ddot{\bf r}_{0}^{T}{\bf R}_{0,q_{i}}{\bf r}_{1}+{\bf r}_{0,q_{i}}^{T}\ddot{\bf R}_{0}{\bf r}_{1}+{\bf r}_{1}^{T}{\bf R}_{0,q_{i}}^{T}\ddot{\bf R}_{0}{\bf r}_{1}\right)d\mathcal{V}}}\\ {{\displaystyle F_{c,i}^{1}=\int_{\mathcal{V}}\rho\left(\ddot{\bf r}_{0}^{T}{\bf R}_{0}{\bf r}_{1,q_{i}}+{\bf r}_{1}^{T}\ddot{\bf R}_{0}^{T}{\bf R}_{0}{\bf r}_{1,q_{i}}\right)d\mathcal{V}}}\end{array} +$$ + +which can be rewritten as + +$$ +\begin{array}{r l}&{F_{c,i}^{0}=\!M\mathbf{r}_{0,q_{i}}^{T}\ddot{\mathbf{r}}_{0}+\left(\ddot{\mathbf{r}}_{0}^{T}\mathbf{R}_{0,q_{i}}+\mathbf{r}_{0,q_{i}}^{T}\ddot{\mathbf{R}}_{0}\right)M\mathbf{r}_{c g}+\left(\mathbf{R}_{0,q_{i}}^{T}\ddot{\mathbf{R}}_{0}\right):\mathbf{l}_{b a s e}}\\ &{F_{c,i}^{1}=\!\ddot{\mathbf{r}}_{0}^{T}\mathbf{R}_{0}M\mathbf{r}_{c g,q_{i}}+\left(\ddot{\mathbf{R}}_{0}^{T}\mathbf{R}_{0}\right):\mathsf{A}_{b a s e,i}}\end{array} +$$ + +using (E.7) and the volume integral notations introduced above. + +# Centrifugal stiffness matrix + +Inserting (1.12) into (1.9), expanding and sorting the terms into the $_{2\times2}$ block matrix form (1.19) yield the centrifugal stiffness matrix elements + +$$ +\begin{array}{l}{{k_{c,i j}^{00}=\displaystyle\int_{\mathcal{V}}\rho\left(\mathbf{r}_{0,q_{i}q_{j}}^{T}\tilde{\mathbf{r}}_{0}+\mathbf{r}_{0,q_{i}}^{T}\tilde{\mathbf{r}}_{0,q_{j}}+\mathbf{r}_{0,q_{i},q_{j}}^{T}\tilde{\mathbf{R}}_{0}\mathbf{r}_{1}+\breve{\mathbf{r}}_{0}^{T}\mathbf{R}_{0,q_{i},q_{j}}\mathbf{r}_{1}+\breve{\mathbf{r}}_{0,q_{j}}^{T}\mathbf{R}_{0,q_{i}}\mathbf{r}_{1}+\mathbf{r}_{0,q_{i}}^{T}\breve{\mathbf{R}}_{0,q_{j}}\mathbf{r}_{1}+\right.}}\\ {{\left.\qquad\quad+\mathbf{\sigma}_{1}^{T}\mathbf{R}_{0,q_{i}}^{T}\tilde{\mathbf{R}}_{0,q_{j}}\mathbf{r}_{1}+\mathbf{r}_{1}^{T}\mathbf{R}_{0,q_{i},q_{j}}^{T}\tilde{\mathbf{R}}_{0}\mathbf{r}_{1}\right)d\mathcal{V}}}\\ {{k_{c,i j}^{01}=\displaystyle\int_{\mathcal{V}}\rho\left(\breve{\mathbf{r}}_{0}^{T}\mathbf{R}_{0,q_{j}}^{T}\mathbf{r}_{1,q_{i}}+\breve{\mathbf{r}}_{0,q_{j}}^{T}\mathbf{R}_{0}\mathbf{r}_{1,q_{i}}+\mathbf{r}_{0,q_{i}}^{T}\breve{\mathbf{R}}_{0}\mathbf{r}_{1,q_{j}}+\breve{\mathbf{r}}_{0}^{T}\mathbf{R}_{0,q_{i}}\mathbf{r}_{1,q_{j}}\right.}}\\ {{\left.\qquad\quad+\mathbf{\sigma}_{1}^{T}\breve{\mathbf{R}}_{0}^{T}\mathbf{R}_{0,q_{j}}\mathbf{r}_{1,q_{i}}+\mathbf{r}_{1}^{T}\breve{\mathbf{R}}_{0,q_{j}}^{T}\mathbf{R}_{0}\mathbf{r}_{1,q_{i}}+\mathbf{r}_{1}^{T}\breve{\mathbf{R}}_{0}^{T}\mathbf{R}_{0,q_{i}}\mathbf{r}_{1,q_{j +$$ + +which can be rewritten as + +$$ +\begin{array}{r l}&{k_{c,i j}^{00}=M\left(\pmb{\Gamma}_{0,q_{i},q_{j}}^{T}\ddot{\mathbf{p}}_{0}+\pmb{\Gamma}_{0,q_{i}}^{T}\ddot{\mathbf{p}}_{0,q_{j}}\right)+\left(\pmb{\Gamma}_{0,q_{i},q_{j}}^{T}\ddot{\mathbf{R}}_{0}+\ddot{\mathbf{r}}_{0}^{T}\pmb{\mathsf{R}}_{0,q_{i},q_{j}}+\ddot{\mathbf{r}}_{0,q_{j}}^{T}\pmb{\mathsf{R}}_{0,q_{i}}+\pmb{\Gamma}_{0,q_{i}}^{T}\ddot{\mathbf{R}}_{0,q_{j}}\right)M\mathbf{r}_{c g}}\\ &{\quad\quad\quad+\left(\pmb{\mathsf{R}}_{0,q_{i}}^{T}\ddot{\mathbf{R}}_{0,q_{j}}+\pmb{\mathsf{R}}_{0,q_{i},q_{j}}^{T}\ddot{\mathbf{R}}_{0}\right):\mathsf{I}_{b a s e}}\\ &{k_{c,i j}^{01}=\left(\ddot{\pmb{\Gamma}}_{0,q_{j}}^{T}\mathbf{R}_{0}+\ddot{\mathbf{r}}_{0}^{T}\mathbf{R}_{0,q_{j}}\right)M\mathbf{r}_{c g,q_{i}}+\left(\pmb{\Gamma}_{0,q_{i}}^{T}\ddot{\mathbf{R}}_{0}+\ddot{\mathbf{r}}_{0}^{T}\mathbf{R}_{0,q_{i}}\right)M\mathbf{r}_{c g,q_{j}}}\\ &{\quad\quad\quad+\left(\ddot{\pmb{\mathsf{R}}}_{0}^{T}\mathbf{R}_{0,q_{j}}+\ddot{\mathbf{R}}_{0,q_{j}}^{T}\mathbf{R}_{0}\right):\mathsf{A}_{b a s e,i}+\left(\ddot{\pmb{\mathsf{R}}}_{0}^{T}\mathbf{R}_{0,q_{i}}+\pmb{\mathsf{R}}_{0,q_{i}}^{T}\ddot{\mathbf{R}}_{0}\right):\mathsf{A}_{b a s e,j}}\\ &{k_{c,i j}^{11}=\ddot{\pmb{\Gamma}}_{0}^{T}\mathbf{R}_{0}M\mathbf{r}_{c g,q_{i},q_{j}}+\left(\pmb{\mathsf{R}}_{0}^{T}\ +$$ + +using (E.7) and the volume integral notations introduced above. + +# Generic implementation of inertia forces from a substructure + +For the implementation of the above inertia force components, each substructure object must provide a function computing all scalers, vectors and matrices marked in blue in Equations (1.22), (1.29), (1.31), (1.33) and (1.35). Some have physical meaning, e.g. the total mass of the substructure $M$ , the current center of gravity $\boldsymbol{\mathsf{r}}_{c g}$ and its DOF derivatives to the first and second order $(\mathsf{r}_{c g,q_{i}}\ a n d\ r_{c g,q_{i}q_{j}})$ , and the current matrix of rotational moments of inertia components $\mid_{b a s e}$ as given by (1.23). The remaining scalars can be reduced to these two unique volume integrals over the substructure that are the entries of the local mass and nonlinear gyroscopic matrices: + +$$ +m_{i j}^{11}=\int_{\mathcal{V}}\rho\,\mathbf{r}_{1,q_{i}}^{T}\mathbf{r}_{1,q_{j}}d\mathcal{V}\;\;\mathbf{a}\mathsf{n}\mathsf{d}\;\;h_{i j k}^{111}=\int_{\mathcal{V}}\rho\,\mathbf{r}_{1,q_{i}}^{T}\mathbf{r}_{1,q_{j}q_{k}}d\mathcal{V} +$$ + +for all $i,j,k\,\in\,{\bf d}_{b}$ where ${\mathfrak{d}}_{b}$ is the DOF index vector for the substructure $b$ . The remaining matrices can be reduced to these three unique volume integrals over the substructure: + +$$ +\mathsf{A}_{b a s e,i}\equiv\int_{\mathcal{V}}\mathsf{r}_{1,q_{i}}\mathsf{r}_{1}^{T}d\mathcal{V}\;,\;\;\mathsf{A}_{b a s e,1,i j}\equiv\int_{\mathcal{V}}\rho\,\mathsf{r}_{1,q_{j}}\mathsf{r}_{1,q_{i}}^{T}d\mathcal{V}\;\;\mathrm{and}\;\;\mathsf{A}_{b a s e,2,i j}\equiv\int_{\mathcal{V}}\rho\,\mathsf{r}_{1,q_{i},q_{j}}\mathsf{r}_{1}^{T}d\mathcal{V} +$$ + +for all $i,j\in{\bf d}_{b}$ . Note that the following symmetry rules apply: ${\pmb{\mathsf{A}}}_{b a s e,2,j i}={\pmb{\mathsf{A}}}_{b a s e,2,i j}$ and $\mathbf{A}_{b a s e,1,j i}=\mathbf{A}_{b a s e,1,i j}^{T}$ + +All scalars, vectors, and matrices for each substructure must be computed in every time step for all combinations of substructure DOFs. The combinations leads to a large number of computations, e.g. the nonlinear gyroscopic elements combines over three DOF indices, with symmetry for two of them, resulting in $\left(N_{b}^{3}+N_{b}^{2}\right)/2$ scalar values (where $N_{b}$ is the number of DOFs in structure $b$ ) to be computed for the nonlinear gyroscopic matrix. The computations of (1.23), (1.36), (1.37), and the derivatives of $\pmb{r}_{c g}$ in the object function of the substructure must therefore be optimized for speed. + +Each substructure object must also provide a function computing its contributions to the positions and orienta­ tions of any connection points to subsequent substructures. + +# 1.1.3 Gravitational forces on and from a substructure + +# 1.1.4 Linear elastic forces inside a substructure + +The flexible substructures will have internal forces that can be modeled as linear elastic forces. These forces can be included in the Lagrange equation through their potential energy $V_{e,b}=V_{e,b}(\mathbf{q}_{b})$ which is only a function of the substructure DOFs. Thus, the internal forces and the associated tangent stiffness matrix do not lead to coupling to the other substructures as the inertia forces. + +# 1.1.5 Linear purely dissipative forces inside a substructure + +The dissipation of vibrational energy in a substructure ... + +# 1.1.6 Generalized forces from external forces and moments + +... aerodynamic forces ... + +# 1.2 Substructure models + +The following model types can be used to model substructures: + +• Rigid body + +• Flexible co­rotational finite beam element body • Flexible linear Craig­Bampton super­element body + +The following sections contain + +# 1.2.1 Rigid body substructure + +The position vector of particles on the rigid body is purely a function of the local spacial coordinates + +$$ +\mathbf{r}_{1,b}=\left\{\begin{array}{l}{x}\\ {y}\\ {z}\end{array}\right\} +$$ + +where $(x,y,z)\in\mathcal{V}$ + +# 1.2.2 Co­rotational beam substructure + +The nonlinear co­rotational finite beam element methodology used for this type of substructures is described in details in Appendix A. The methodology is similar to Crisfield [2] and Krenk [3] except that the nonlinear geometric formulation is explicit, i.e., the element coordinate systems and the local small rotations of the nodes are given as explicit nonlinear functions of the global nodal DOFs. The elastic deformation and compliance of the element is adapted from the equilibrium element proposed by Krenk and Couturier [4]. + +A co­rotational beam substructure consists of finite beam elements with two end­nodes having six DOFs: three translations and three rotations (using Rodrigues parameters) in the substructure coordinate system. The DOF vector of the substructure + +# Inertia forces + +The volume integral over a co­rotational beam substructure can be computed as the sum of integrals over the volume $\mathcal{V}_{n}$ of element $n$ , here generically written as + +$$ +\int_{\mathcal{V}}\left(\mathbf{\Phi}\right)d\mathcal{V}=\sum_{n=1}^{N_{e,b}}\int_{\mathcal{V}_{n}}\left(\mathbf{\Phi}\right)d\mathcal{V}_{n}=\sum_{n=1}^{N_{e,b}}\left(\frac{l_{n}}{2}\int_{-1}^{1}\int_{\mathcal{A}}\left(\mathbf{\Phi}\right)d\mathcal{A}d\zeta\right) +$$ + +where each volume integral over the element of initial length $l_{n}$ are split into an area integral over each cross­ section and a line integral over the non­dimensional element coordinate $\zeta$ from the mid­point to the end nodes at $\zeta=\pm1$ . + +The position vector of particles on element number $n$ is + +$$ +\mathbf{r}_{1,b,n}=\mathbf{r}_{\mathrm{mid},b,n}\left(\mathbf{q}_{b,n}\right)+\mathbf{E}_{b,n}\left(\mathbf{q}_{b,n}\right)\mathbf{v}_{b,n}\left(\mathbf{q}_{b,n};x,y,\zeta\right) +$$ + +where ${\boldsymbol{\mathsf{r}}}_{{\boldsymbol{\mathsf{m i d}}},{\boldsymbol{b}},n}$ is a vector from the substructure base to the mid­point of the element linearly dependent on the displacement DOFs of element nodes, $\mathsf{E}_{b,n}$ is the element coordinate system dependent on all twelve nodal DOFs $\mathbf{q}_{b,n}$ of element $n$ , and the cross­sectional displacement vector is given by + +$$ +\begin{array}{r}{\mathbf{v}_{b,n}=\left\{\begin{array}{c}{x}\\ {y}\\ {\frac{1}{2}l_{n}\zeta}\end{array}\right\}+\left\{\begin{array}{c}{u_{x,n}(\zeta)}\\ {u_{y,n}(\zeta)}\\ {u_{z,n}(\zeta)}\end{array}\right\}+\left[\begin{array}{c c c}{0}&{-\theta_{z,n}(\zeta)}&{\theta_{y,n}(\zeta)}\\ {\theta_{z,n}(\zeta)}&{0}&{-\theta_{x,n}(\zeta)}\\ {-\theta_{y,n}(\zeta)}&{\theta_{x,n}(\zeta)}&{0}\end{array}\right]\left\{\begin{array}{c}{x}\\ {y}\\ {0}\end{array}\right\}}\end{array} +$$ + +where $x,y$ are the coordinates of the cross­section and the cross­sectional translations and (small) rotations are given by the shape function polynomials + +$$ +\left\{\begin{array}{l}{u_{x,n}}\\ {u_{y,n}}\\ {u_{z,n}}\\ {\theta_{x,n}}\\ {\theta_{y,n}}\\ {\theta_{z,n}}\end{array}\right\}=\sum_{p=0}^{P+3}\mathbf{N}_{n,p}\,\mathbf{g}\left(\mathbf{q}_{b,n}\right)\zeta^{p}=\sum_{p=0}^{P+3}\left[\begin{array}{l}{\bar{\mathbf{N}}_{n,p}}\\ {\tilde{\mathbf{N}}_{n,p}}\end{array}\right]\mathbf{g}\left(\mathbf{q}_{b,n}\right)\zeta^{p} +$$ + +where $\aleph_{n,p}$ are the $6\uptimes7$ coefficient matrices of the local shape functions, $P$ is the polynomial order of the element properties (e.g. for prismatic elements $P=0$ ), and $\mathfrak{g}$ is the nonlinear $7\times1$ vector function of the twelve nodal DOFs. In the following, we split the local shape function coefficient matrix into two, one $3\mathrm{x7}$ matrix for the displacements $\bar{\aleph}_{n,p}$ and another 3x7 matrix for the rotations $\tilde{\mathsf{N}}_{n,p}$ . + +This form of the cross­sectional displacement vector is inconvenient for the isolation of the spacial variables $x,y,\zeta$ . We therefore rewrite it as + +$$ +\begin{array}{r}{\mathbf{v}_{b,n}=\left\{\begin{array}{c}{x}\\ {y}\\ {\frac{1}{2}l_{n}\zeta}\end{array}\right\}+\displaystyle\sum_{p=0}^{P+3}\left\{\begin{array}{c}{u_{x,n,p}}\\ {u_{y,n,p}}\\ {u_{z,n,p}}\end{array}\right\}\zeta^{p}+\displaystyle\sum_{p=0}^{P+3}\left(\left\{\begin{array}{c}{0}\\ {\theta_{z,n,p}}\\ {-\theta_{y,n,p}}\end{array}\right\}x+\left\{\begin{array}{c}{-\theta_{z,n,p}}\\ {0}\\ {\theta_{x,n,p}}\end{array}\right\}y\right)\zeta^{p}}\end{array} +$$ + +where the polynomial coefficients of the cross­sectional displacements and rotations are given by + +$$ +\left\{\begin{array}{c}{u_{x,n,p}}\\ {u_{y,n,p}}\\ {u_{z,n,p}}\end{array}\right\}=\bar{\aleph}_{n,p}\,\mathbf{g}\left(\mathbf{q}_{b,n}\right)\quad\mathsf{a n d}\quad\left\{\begin{array}{c}{\theta_{x,n,p}}\\ {\theta_{y,n,p}}\\ {\theta_{z,n,p}}\end{array}\right\}=\tilde{\aleph}_{n,p}\mathbf{g}\left(\mathbf{q}_{b,n}\right) +$$ + +that depend on the nodal DOFs of the element. + +Using (1.43) with (1.42), the position vector for element $n$ on substructure $b$ (1.40) can be rewritten as + +$$ +\mathbf{r}_{1,b,n}=\sum_{p=0}^{P+3}\left(\mathbf{r}_{o,b,n,p}+x\,\mathbf{r}_{x,b,n,p}+y\,\mathbf{r}_{y,b,n,p}\right)\zeta^{p} +$$ + +where the sub­vectors are + +$$ +\begin{array}{r l}&{\mathbf{r}_{o,b,n,p}\left(\mathbf{q}_{b,n}\right)=\pmb{\Sigma}_{b,n}\left(\mathbf{q}_{b,n}\right)\,\bar{\mathbf{N}}_{n,p}\,\mathbf{g}\left(\mathbf{q}_{b,n}\right)+\delta_{0p}\,\mathbf{r}_{\mathrm{mid},b,n}\left(\mathbf{q}_{b,n}\right)+\delta_{1p}\,\frac{l_{n}}{2}\mathbf{e}_{3,b,n}\left(\mathbf{q}_{b,n}\right)}\\ &{\mathbf{r}_{x,b,n,p}\left(\mathbf{q}_{b,n}\right)=\pmb{\Sigma}_{b,n}\left(\mathbf{q}_{b,n}\right)\,\mathbf{P}_{x}\,\tilde{\mathbf{N}}_{n,p}\,\mathbf{g}\left(\mathbf{q}_{b,n}\right)+\delta_{0p}\,\mathbf{e}_{1,b,n}\left(\mathbf{q}_{b,n}\right)}\\ &{\mathbf{r}_{y,b,n,p}\left(\mathbf{q}_{b,n}\right)=\pmb{\Sigma}_{b,n}\left(\mathbf{q}_{b,n}\right)\,\mathbf{P}_{y}\,\tilde{\mathbf{N}}_{n,p}\,\mathbf{g}\left(\mathbf{q}_{b,n}\right)+\delta_{0p}\,\mathbf{e}_{2,b,n}\left(\mathbf{q}_{b,n}\right)}\end{array} +$$ + +where $\delta_{i j}$ is Kronecker’s delta which is 1 for $i=j$ , and otherwise 0. Matrices $\mathsf{P}_{x}$ and $\mathsf{P}_{y}$ are constant permutation matrices: + +$$ +\begin{array}{r}{\mathsf{P}_{x}=\left[\begin{array}{c c c}{0}&{0}&{0}\\ {0}&{0}&{1}\\ {0}&{-1}&{0}\end{array}\right]\quad\mathsf{a n d}\quad\mathsf{P}_{y}=\left[\begin{array}{c c c}{0}&{0}&{-1}\\ {0}&{0}&{0}\\ {1}&{0}&{0}\end{array}\right]}\end{array} +$$ + +Note that $\mathsf{P}_{x}$ and $\mathsf{P}_{y}$ are equal to the matrix $\boldsymbol{\mathsf{N}}$ in Eq. (E.3) related to the angular derivative of a rotation matrix with a constant unit­vectors in the $x$ ­ and $y$ ­directions, respectively. + +In the following, the element­index $n$ and the substructure­index $b$ are omitted for brevity. + +Integration over the cross­sectional area spanned by $x$ and $y$ of the element coordinate system are defined by the mass per unit­length $m$ , the center of gravity coordinates $x_{c g}$ and $y_{c g}$ , and the moments of inertia about the $x$ ­ and $y$ ­axes and their cross­coupling. These cross­sectional properties are given by polynomials to the order of $P$ along the element lengthwise coordinate $\zeta\in[-1,1]$ as + +$$ +\begin{array}{l}{\displaystyle\int_{A}\rho\,d A=\sum_{r=0}^{P}a_{m,r}\,\zeta^{r}\;,\;\;\int_{A}\rho x\,d A=\sum_{r=0}^{P}a_{m x_{c g},r}\,\zeta^{r}\;,\;\;\int_{A}\rho y\,d A=\sum_{r=0}^{P}a_{m y_{c g},r}\,\zeta^{r}\;,}\\ {\displaystyle\int_{A}\rho x^{2}\,d A=\sum_{r=0}^{P}a_{I_{x x},r}\,\zeta^{r}\;,\;\;\displaystyle\int_{A}\rho y^{2}\,d A=\sum_{r=0}^{P}a_{I_{y y},r}\,\zeta^{r}\;,\;\;\displaystyle\int_{A}\rho x y\,d A=\sum_{r=0}^{P}a_{I_{x y},r}\,\zeta^{r}}\end{array} +$$ + +where the polynomial coefficients are generated from the model input in a pre­simulation processing step. The coefficients $a_{m x_{c g},r}$ and $a_{m y_{c g},r}$ represent a polynomial fit to the product of the mass per unit­length $m$ and the individual center of gravity coordinates $x_{c g}$ and $y_{c g}$ in the element coordinate system. + +Integration over an element with the local position vector (1.45) and using (1.48), the mass of the element can be computed as + +$$ +M=\frac{l}{2}\sum_{r=0}^{P}c(r)a_{m,r} +$$ + +where the coefficient function is given by + +$$ +c\left(r\right)={\frac{\left(-1\right)^{r}+1}{1+r}}={\left\{\begin{array}{l l}{2/(1+r)}&{r\;\;{\mathsf{e v e n}}}\\ {\;\;\;\;0}&{r\;\;\;{\mathsf{o d d}}}\end{array}\right.} +$$ + +The element center of gravity position times the element mass and the nodal derivatives can be computed as + +$$ +\begin{array}{l l}{{\displaystyle M\,{\bf r}_{c g}=\frac{l}{2}\sum_{p=0}^{P+3}\left(\sum_{r=0}^{P}c(p+r)\left(a_{m,r}\,{\bf r}_{o,p}+a_{m x_{c g},r}\,{\bf r}_{x,p}+a_{m y_{c g},r}\,{\bf r}_{y,p}\right)\right)}}\\ {{\displaystyle M\,{\bf r}_{c g,q_{i}}=\frac{l}{2}\sum_{p=0}^{P+3}\left(\sum_{r=0}^{P}c(p+r)\left(a_{m,r}\,{\bf r}_{o,p,q_{i}}+a_{m x_{c g},r}\,{\bf r}_{x,p,q_{i}}+a_{m y_{c g},r}\,{\bf r}_{y,p,q_{i}}\right)\right)}}\\ {{\displaystyle M\,{\bf r}_{c g,q_{i},q_{j}}=\frac{l}{2}\sum_{p=0}^{P+3}\left(\sum_{r=0}^{P}c(p+r)\left(a_{m,r}\,{\bf r}_{o,p,q_{i},q_{j}}+a_{m x_{c g},r}\,{\bf r}_{x,p,q_{i},q_{j}}+a_{m y_{c g},r}\,{\bf r}_{y,p,q_{i},q_{j}}\right)\right)}}\end{array} +$$ + +where $c(p+r)$ is given by (1.49). The scalars and matrices of (1.36) and (1.37) have the generic form + +$$ +\frac{l}{2}\sum_{p=0}^{P+3}\left(\sum_{r=0}^{P}\left(\sum_{q=0}^{P+3}c\left(q+r+p\right)\mathcal{G}\left\{\mathbf{a}_{r},\mathbf{u}_{o,q},\mathbf{u}_{x,q},\mathbf{u}_{y,q},\mathbf{w}_{o,p},\mathbf{w}_{x,p},\mathbf{w}_{y,p}\right\}\right)\right) +$$ + +where the generic operator is defined as + +$$ +\begin{array}{r}{\mathcal{G}\left\{\mathbf{a},\mathbf{u}_{o},\mathbf{u}_{x},\mathbf{u}_{y},\mathbf{w}_{o},\mathbf{w}_{x},\mathbf{w}_{y}\right\}={a}_{m}\,\mathbf{u}_{o}^{T}\mathbf{w}_{o}+{a}_{m x_{c g}}\,\left(\mathbf{u}_{o}^{T}\mathbf{w}_{x}+\mathbf{u}_{x}^{T}\mathbf{w}_{o}\right)+{a}_{m y_{c g}}\,\left(\mathbf{u}_{o}^{T}\mathbf{w}_{y}+\mathbf{u}_{y}^{T}\mathbf{w}_{o}\right)}\\ {+\;{a}_{i_{x x}}\,\mathbf{u}_{x}^{T}\mathbf{w}_{x}+{a}_{i_{x y}}\,\left(\mathbf{u}_{x}^{T}\mathbf{w}_{y}+\mathbf{u}_{y}^{T}\mathbf{w}_{x}\right)+{a}_{i_{y y}}\,\mathbf{u}_{y}^{T}\mathbf{w}_{y}\,\,\,}\end{array} +$$ + +where $\mathbf{a}=[a_{m},a_{m x_{c g}},a_{m y_{c g}},a_{i_{x x}},a_{i_{x y}},a_{i_{y y}}]$ is a list with the polynomial coefficients of the material properties (such that ${\bf a}_{r}$ contains the $r$ ’th order coefficients), and the vectors $\mathbf{u}_{o}$ , $\mathbf{u}_{x}$ , $\mathbf{u}_{y}$ , $\pmb{{\mathsf{w}}}_{o}$ , $\mathbf{w}_{x}$ , and $\pmb{w}_{y}$ are 3x1 column + +vectors or $\mathsf{1x3}$ row vectors for scalar or matrix evaluations, respectively. The local mass matrix contribution from an element (the entries of the element mass matrix) is + +$$ +m_{i j}^{11}=\frac{l}{2}\sum_{p=0}^{P+3}\left(\sum_{r=0}^{P}\left(\sum_{q=0}^{P+3}c\left(q+r+p\right)\;{\mathcal G}\left\{\mathbf a_{r},\mathbf r_{o,q,q_{i}},\mathbf r_{x,q,q_{i}},\mathbf r_{y,q,q_{i}},\mathbf r_{o,p,q_{j}},\mathbf r_{x,p,q_{j}},\mathbf r_{y,p,q_{j}}\right\}\right)\right) +$$ + +and its contribution to the local nonlinear gyroscopic matrix (the nonlinear element gyroscopic matrix) is + +$$ +h_{i j k}^{111}=\frac{l}{2}\sum_{p=0}^{P+3}\left(\sum_{r=0}^{P}\left(\sum_{q=0}^{P+3}c\left(q+r+p\right)\;\mathcal{G}\left\{\mathbf{a}_{r},\mathbf{r}_{o,q,q_{i}},\mathbf{r}_{x,q,q_{i}},\mathbf{r}_{y,q,q_{i}},\mathbf{r}_{o,p,q_{j},q_{k}},\mathbf{r}_{x,p,q,q_{j},q_{k}},\mathbf{r}_{y,p,q,q_{k}}\right\}\right)\right) +$$ + +The contribution of an element to the matrix related to the rotational inertia of the substructure about its base is + +$$ +\mathbf{l}_{b a s e}=\frac{l}{2}\sum_{p=0}^{P+3}\left(\sum_{r=0}^{P}\left(\sum_{q=0}^{P+3}c\left(q+r+p\right)\,\mathcal{G}\left\{\mathbf{a}_{r},\mathbf{r}_{o,q}^{T},\mathbf{r}_{x,q}^{T},\mathbf{r}_{y,q}^{T},\mathbf{r}_{o,p}^{T},\mathbf{r}_{x,p}^{T},\mathbf{r}_{y,p}^{T}\right\}\right)\right) +$$ + +Finally, the matrices needed for the computation of the remaining inertia forces from an element are + +$$ +\begin{array}{r l}&{\displaystyle\mathsf{A}_{b a s e,i,i}=\frac{l}{2}\sum_{p=0}^{P+3}\left(\sum_{r=0}^{P}\left(\sum_{q=0}^{P+3}c(q+r+p)\;\mathcal{G}\left\{\mathbf{a}_{r},\mathbf{r}_{o,p,q_{i}}^{T},\mathbf{r}_{x,p,q_{i}}^{T},\mathbf{r}_{y,p,q_{i}}^{T},\mathbf{r}_{o,q}^{T},\mathbf{r}_{x,q}^{T},\mathbf{r}_{y,q}^{T}\right\}\right)\right)}\\ &{\displaystyle\mathsf{A}_{b a s e,1,i,j}=\frac{l}{2}\sum_{p=0}^{P+3}\left(\sum_{r=0}^{P}\left(\sum_{q=0}^{P+3}c(q+r+p)\;\mathcal{G}\left\{\mathbf{a}_{r},\mathbf{r}_{o,p,q_{j}}^{T},\mathbf{r}_{x,p,q_{j}}^{T},\mathbf{r}_{y,p,q_{j}}^{T},\mathbf{r}_{o,q_{i}}^{T},\mathbf{r}_{x,q,q_{i}}^{T},\mathbf{r}_{y,q,\bar{q}_{i}}^{T}\right\}\right)\right)}\\ &{\displaystyle\mathsf{A}_{b a s e,2,i,j}=\frac{l}{2}\sum_{p=0}^{P+3}\left(\sum_{r=0}^{P}\left(\sum_{q=0}^{P+3}c(q+r+p)\;\mathcal{G}\left\{\mathbf{a}_{r},\mathbf{r}_{o,p,q_{i},q_{j}}^{T},\mathbf{r}_{x,p,q_{i},q_{j}}^{T},\mathbf{r}_{y,p,q_{i},q_{j}}^{T},\mathbf{r}_{x,q}^{T},\mathbf{r}_{x,q}^{T},\mathbf{r}_{y,q}^{T},\mathbf{r}_{y,q}^{T}\right\}\right)\right)}\end{array} +$$ + +The generic form of these contributions to the inertia forces from a single element can be implemented in a single object function that only need the polynomial coefficient sub­vectors (1.46) and their first nodal derivatives + +$$ +\begin{array}{r l}&{\mathbf{r}_{o,b,n,p,q_{i}}=\!\!\mathsf{E}_{b,n,q_{i}}\,\bar{\mathbf{N}}_{n,p}\,\mathbf{g}+\mathsf{E}_{b,n}\,\bar{\mathbf{N}}_{n,p}\,\mathbf{g}_{q_{i}}+\delta_{0p}\,\mathbf{r}_{\mathrm{mid},b,n,q_{i}}+\delta_{1p}\,\frac{l_{n}}{2}\mathbf{e}_{3,b,n,q_{i}}}\\ &{\mathbf{r}_{x,b,n,p,q_{i}}=\!\!\mathsf{E}_{b,n,q_{i}}\,\mathsf{P}_{x}\,\tilde{\mathbf{N}}_{n,p}\,\mathbf{g}+\mathsf{E}_{b,n}\,\mathsf{P}_{x}\,\tilde{\mathbf{N}}_{n,p}\,\mathbf{g}_{,q_{i}}+\delta_{0p}\,\mathbf{e}_{1,b,n,q_{i}}}\\ &{\mathbf{r}_{y,b,n,p,q_{i}}=\!\!\mathsf{E}_{b,n,q_{i}}\,\mathsf{P}_{y}\,\tilde{\mathbf{N}}_{n,p}\,\mathbf{g}+\mathsf{E}_{b,n}\,\mathsf{P}_{y}\,\tilde{\mathbf{N}}_{n,p}\,\mathbf{g}_{,q_{i}}+\delta_{0p}\,\mathbf{e}_{2,b,n,q_{i}}}\end{array} +$$ + +and second nodal derivatives + +$$ +\begin{array}{r l}{\mathbf{r}_{o,b,n,p,q_{i},q_{j}}=\mathbf{E}_{b,n,q_{i},q_{j}}\,\bar{\mathbf{N}}_{n,p}\,\mathbf{g}+\mathbf{E}_{b,n,q_{j}}\,\bar{\mathbf{N}}_{n,p}\,\mathbf{g}_{q_{i}}}&{}\\ {+\,\mathbf{E}_{b,n,q_{i}}\,\bar{\mathbf{N}}_{n,p}\,\mathbf{g}_{q_{j}}+\mathbf{E}_{b,n}\,\bar{\mathbf{N}}_{n,p}\,\mathbf{g}_{q_{i},q_{j}}+\delta_{1p}\,\frac{l_{n}}{2}\mathbf{e}_{3,b,n,q_{i},q_{j}}}&{}\\ {\mathbf{r}_{x,b,n,p,q_{i},q_{j}}=\mathbf{E}_{b,n,q_{i},q_{j}}\,\mathbf{p}_{x}\,\tilde{\mathbf{N}}_{n,p}\,\mathbf{g}+\mathbf{E}_{b,n,q_{j}}\,\mathbf{p}_{x}\,\tilde{\mathbf{N}}_{n,p}\,\mathbf{g}_{q_{i}}}&{}\\ {+\,\mathbf{E}_{b,n,q_{i}}\,\mathbf{p}_{x}\,\tilde{\mathbf{N}}_{n,p}\,\mathbf{g}_{q_{j}}+\mathbf{E}_{b,n}\,\mathbf{p}_{x}\,\tilde{\mathbf{N}}_{n,p}\,\mathbf{g}_{q_{i},q_{j}}+\delta_{0p}\,\mathbf{e}_{1,b,n,q_{i},q_{j}}}&{}\\ {\mathbf{r}_{y,b,n,p,q_{i},q_{j}}=\mathbf{E}_{b,n,q_{i},q_{j}}\,\mathbf{p}_{y}\,\tilde{\mathbf{N}}_{n,p}\,\mathbf{g}+\mathbf{E}_{b,n,q_{j}}\,\mathbf{p}_{y}\,\tilde{\mathbf{N}}_{n,p}\,\mathbf{g}_{q_{i}}}&{}\\ {+\,\mathbf{E}_{b,n,q_{i}}\,\mathbf{p}_{y}\,\tilde{\mathbf{N}}_{n,p}\,\mathbf{g}_{q_{j}}+\mathbf{E}_{b,n}\,\mathbf{p}_{y}\,\tilde{\mathbf{N}}_{n,p}\,\mathbf{g}_{q_{i},q_{j}}+\delta_{0p}\,\mathbf{e}_{2,b,n,q_{i},q_ +$$ + +where the element­index $n$ and substructure­index $b$ are included again. Note that the coefficient sub­vectors (1.46) and its first and second derivatives (1.57) and (1.58) must be evaluated at the current nodal DOFs $\mathbf{q}_{b,n}$ of the element. The $\mathsf{3}\mathsf{x}\mathsf{3}$ element coordinate system $\mathsf{E}_{b,n}$ , the $7\times1$ vector function g, and their derivatives are derived as explicit functions of the nodal DOFs in Appendix A. The mid­element position vector is simply a linear function of the $_{2\times3}$ displacements of the end­nodes (A.2). + +The substructure object must return the following scalar, vectors, and matrices: + +$M$ Total mass of the substructure which is constant. rcg = nN=e,1b rcg,n The $3\!\times\!1$ position vector for the center of gravity of the substructure where the contributions from each element are given by (1.51a). $\begin{array}{r}{\mathbf{r}_{c g,q_{i}}=\sum_{n=1}^{N_{e,b}}\mathbf{r}_{c g,n,q_{i}}}\end{array}$ The first nodal derivatives of the 3x1 position vector for the center of gravity of the substructure where the contributions from each element are given by (1.51b). Maximum two elements will contribute to each DOF derivative. $\begin{array}{r}{\mathbf{r}_{c g,q_{i},q_{j}}=\sum_{n=1}^{N_{e,b}}\mathbf{r}_{c g,n,q_{i},q_{j}}}\end{array}$ The second nodal derivatives of the 3x1 position vector for the center of gravity of the substructure where the contributions from each element are given by (1.51c). Maximum two elements will contribute to each DOF derivative. $\mathbb{M}^{11}$ The local mass matrix of the co­rotational substructure built up by collecting all 12x12 element mass matrices (1.53) on its diagonal. $\mathsf{H}^{111}$ The local nonlinear gyroscopic matrix of the co­rotational substructure built up by collecting all 12x12x12 element nonlinear gyroscopic matrices (1.54) on its diagonal. This local matrix is three­dimensional where the twelve nodal velocities at the ends of each element cause gyroscopic forces on the twelve DOF equations of the same element. $\begin{array}{r}{\mathbf{l}_{b a s e}=\sum_{n=1}^{N_{e,b}}\mathbf{l}_{b a s e,n}}\end{array}$ The $\mathsf{3}\mathsf{x}\mathsf{3}$ rotational inertia matrix where the contributions from each element are given by (1.55). + +# Gravity forces + +# Elastic stiffness forces + +Structural damping forces + +# 1.2.3 Linear super­element substructure + +# 2 Rotor aerodynamics + +First in this chapter, a generic formulation is derived for the generalized aerodynamic forces from a blade onto the DOFs of the individual substructures of the structural model. It is important to note that the aerodynamic blade can consist of several substructures, e.g. if the blade has a pitch bearing half way down the blade, also called partial pitch. Second, the aerodynamic forces in a blade cross­section are derived from a generic formulation of the local relative flow in the velocity triangles and the explicit reformulation of the unsteady airfoil aerodynamic model developed in [5]. Finally, the aerodynamic sub­model of the far wake induction used in CASEStab is described based on the dynamic inflow model with localized induction developed for HAWC2 [6]. + +# 2.1 Generalized aerodynamic forces on and from a blade + +The contribution to the non­conservative generalized force on coordinate $q_{i}$ from the aerodynamic pressure P on the surface $\boldsymbol{S}$ of blade number $\beta$ (not to be confused with the substructure index $b$ ) can be written as + +$$ +Q_{i,\beta}=\int_{S}\frac{\partial\mathbf{r}_{p,\beta}^{T}}{\partial q_{i}}\mathbf{p}_{\beta}\;d S +$$ + +where $\mathsf{r}_{p,\beta}$ is a position vector of the partial force vector $\mathsf{P}_{\beta}d S$ in the inertial frame. The aerodynamic surface pressure vector $\mathsf{P}_{\beta}$ depend on time $t$ , the displacements $\mathfrak{q}$ and velocities $\dot{\mathbf{q}}$ of the turbine structure (added mass effects are neglected in CASEStab), and the deterministic and stochastic input ${\bf w}(t)$ from the wind field. In the unsteady Blade Element Momentum (BEM) formulation, the aerodynamic forces in the cross­section are described by the predefined curve integrations of the surface pressure around the cross­section along each blade, and the unsteady aerodynamic effects of shed and trailed vorticity, dynamic stall effects, and far wake effects on the induced velocities are handled by different engineering models with aerodynamic state­variables collected in the vector ${\pmb x}_{a}$ (cf. Sections 2.2.2 and 2.4). + +With reference to Figure 2.1, let $\mathbf{v}_{c}=\left\{x,y,0\right\}^{T}$ be the local position vector of the surface of the cross­section described in the chord coordinate system (CCS), where the first axis has the direction of the chord (positive to­ wards leading edge), the second axis is perpendicular to the chord (positive towards suction side), and the third axis is the normal to the plane of the airfoil for which the surface pressure (later the aerodynamic coefficients) is defined. The integration over the entire blade (2.1) can then be rewritten as + +$$ +Q_{i,\beta}\!=\!\int_{0}^{L}\!\!\oint_{C}\frac{\partial\mathbf{r}_{p,\beta}\left(x,y,s,t,\mathbf{q}\right)^{T}}{\partial q_{i}}\mathbf{P}_{\beta}\left(x,y,s,t,\mathbf{q},\dot{\mathbf{q}},\mathbf{x}_{a},\mathbf{w}\right)d x\,d y\,d s +$$ + +where $C$ is the curve around the cross­section, $L$ is the total curve­length of the blade, and $s$ is the spanwise curve coordinate. The position vector of the (non­deformable) cross­section can then be written as + +$$ +\mathbf{r}_{p,\beta}\left(x,y,s,t,\mathbf{q}\right)=\mathbf{r}_{a c,\beta}\left(s,t,\mathbf{q}\right)+\mathbf{E}_{c,\beta}\left(s,t,\mathbf{q}\right)\mathbf{v}_{c} +$$ + +where $\mathfrak{r}_{a c,\beta}$ is the position vector of the aerodynamic center and $\mathsf{E}_{c,\beta}$ is the orientation matrix of the CCS in the inertia frame, which are both functions of the spanwise coordinate, time and structural displacements q. The orientation of the CCS is defined such that the blade surface normal vector lies in the plane spanned by $\mathfrak{e}_{1,c,\beta}$ and $\mathfrak{e}_{2,c,\beta}$ . Thus, the aerodynamic surface pressure vector can be written as + +$$ +\mathsf{P}_{\beta}=p_{\beta}(x,y,s,t,\mathbf{q},\dot{\mathbf{q}},\mathbf{x}_{a},\mathbf{w})\,\mathsf{E}_{c,\beta}\,(s,t,\mathbf{q})\,\mathsf{n}_{c}\,(s,x,y) +$$ + +where $p_{\beta}$ is the scalar value of the surface pressure around the cross­section, and $\mathbf{n}_{c}\:=\:\left\{n_{x,c},n_{y,c},0\right\}^{T}$ is the surface normal vector described in the CCS where the third coordinate is zero. Using (2.3) and (2.4), the + +![](images/ac7e6ec1853ac26e03566d17b33a6f52f65da377c47d049e52eeb6ca59a149a3.jpg) + +Figure 2.1: Schematics of a blade cross­section with the CCS $\mathsf{E}_{c,\beta}$ with its origin at the aerodynamic center in the quarter chord point. The orientation of this coordinate system is defined such that the surface normal vector lies in the plane spanned by first two unit­vectors $\mathfrak{e}_{1,c,\beta}$ and $\mathfrak{e}_{1,c,\beta}$ . + +generalized aerodynamic force (2.2) can be written as + +$$ +Q_{i,\beta}=\int_{0}^{L}\oint_{C}p_{\beta}\left({\frac{\partial\mathbf{r}_{a c,\beta}^{T}}{\partial q_{i}}}\mathbf{E}_{c,\beta}\mathbf{n}_{c}+\mathbf{v}_{c}^{T}{\frac{\partial\mathbf{E}_{c,\beta}^{T}}{\partial q_{i}}}\mathbf{E}_{c,\beta}\mathbf{n}_{c}\right)\,d x\,d y\,d s +$$ + +where the first term relates to the translation of the aerodynamic center and the second term relates to the rotation of the CCS due to the generalized structural coordinate $q_{i}$ . + +Looking at the inner curve integral for the first term, we note that it can be rewritten as + +$$ +\frac{\partial\mathbf{r}_{a c,\beta}^{T}}{\partial q_{i}}\pmb{\mathsf{E}}_{c,\beta}\oint_{C}p_{\beta}\mathbf{n}_{c}\;d x\,d y=\frac{\partial\mathbf{r}_{a c,\beta}^{T}}{\partial q_{i}}\pmb{\mathsf{E}}_{c,\beta}\left\{\begin{array}{c}{f_{x,\beta}}\\ {f_{y,\beta}}\\ {0}\end{array}\right\} +$$ + +where the local aerodynamic forces in the CCS are defined as + +$$ +f_{x,\beta}=\oint_{C}p_{\beta}n_{x,c}\;d x\,d y\quad{\sf a n d}\quad f_{y,\beta}=\oint_{C}p_{\beta}n_{y,c}\;d x\,d y +$$ + +which are functions of lift and drag forces and the geometric angle of attack, as described in Section 2.2. + +Looking at the inner curve integral for the second term of (2.5), we note that the product of the transposed derivative of the orientation matrix with the orientation matrix itself is a skew­symmetric matrix + +$$ +\frac{\partial\pmb{\mathsf{E}}_{c,\beta}^{T}}{\partial q_{i}}\pmb{\mathsf{E}}_{c,\beta}=\left[\begin{array}{c c c}{0}&{\omega_{3,\beta,i}}&{-\omega_{2,\beta,i}}\\ {-\omega_{3,\beta,i}}&{0}&{\omega_{1,\beta,i}}\\ {\omega_{2,\beta,i}}&{-\omega_{1,\beta,i}}&{0}\end{array}\right] +$$ + +where $\omega_{j,\beta,i}\ (j=1,2,3)$ are the components of the pseudo vector that describes the local rotation of the CCS due to the generalized coordinate $q_{i}$ . Using this property, the inner curve integral for the second term of (2.5) can be written as + +$$ +\oint_{C}p_{\beta}\mathbf{v}_{c}^{T}\frac{\partial\mathbf{E}_{c,\beta}^{T}}{\partial q_{i}}\mathbf{E}_{c,\beta}\mathbf{n}_{c}\;d x\,d y=\oint_{C}p_{\beta}\omega_{3,\beta,i}\left(x\,n_{y,c}-y\,n_{x,c}\right)\,d x\,d y=\omega_{3,\beta,i}M_{\beta} +$$ + +where $M_{\beta}$ is recognized as the aerodynamic moment about $\mathbf{e}_{c,\beta,3}$ in the aerodynamic center defined as + +$$ +M_{\beta}=\oint_{C}p_{\beta}\left(x\,n_{y,c}-y\,n_{x,c}\right)\,d x\,d y +$$ + +Note that only the third component of the pseudo rotation vector $\omega_{\beta,i,3}=\partial\mathbf{e}_{1,c,\beta}^{T}/\partial q_{i}\,\mathbf{e}_{2,c,\beta}$ from variation of the structural coordinate $q_{i}$ is needed, as it describes the rotation of the cross­section about the same axis as the aerodynamic moment. + +Hence, the generalized force on the structural coordinate $q_{i}$ due to aerodynamic forces and moment distribu­ tions on the blade number $\beta$ can be computed from integration along its span as + +$$ +Q_{i,\beta}=\int_{0}^{L}\left(\frac{\partial\mathbf{r}_{a c,\beta}^{T}}{\partial q_{i}}\;\pmb{\mathsf{E}}_{c,\beta}\mathbf{\hat{f}}_{c,\beta}+\frac{\partial\mathbf{e}_{1,c,\beta}^{T}}{\partial q_{i}}\;\pmb{\mathsf{e}}_{2,c,\beta}\;M_{\beta}\right)d s +$$ + +where the local aerodynamic force vector $\mathbf{f}_{c,\beta}\:=\:\{f_{x,\beta},f_{y,\beta},0\}^{T}$ and moment in the blade cross­section are derived from velocity triangles described in the following section. + +# 2.1.1 Generalized aerodynamic forces on and from a substructure of the blade + +The aerodynamic blade number $\beta$ may consist of $N_{\beta}$ substructures with the numbers $b=b_{\beta},\ldots,b_{\beta}+B_{\beta}-$ 1, where $b_{\beta}$ is the number of the first substructure of the blade. The total contribution from blade $\beta$ to the generalized forces is therefore computed as the sum + +$$ +Q_{i,\beta}=\sum_{b=b_{\beta}}^{b_{\beta}+B_{\beta}-1}Q_{i,\beta,b} +$$ + +where $Q_{i,\beta,b}$ is the generalized aerodynamic force on and from the substructure $b$ of the blade $\beta$ . To derive generic expressions for this force, we will write the position of aerodynamic center in a form similar to the structural kinematics Eq. (1.12): + +$$ +\begin{array}{r}{\mathbf{r}_{a c,\beta}\left(t,\mathbf{q};x,y,s\right)=\mathbf{r}_{0,b}\left(t,\mathbf{q}\right)+\mathbf{R}_{0,b}\left(t,\mathbf{q}\right)\mathbf{r}_{a c,1,\beta,b}\left(\mathbf{q}_{b};x,y,s\right)}\end{array} +$$ + +where the translations $\boldsymbol{\mathsf{r}}_{0,b}$ and rotations $\pmb{\mathsf{R}}_{0,b}$ of the base of the substructure $b$ are given by (1.13), and $\mathbf{r}_{a c,1,\beta,b}$ is the local position vector of the aerodynamic center of blade $\beta$ in the frame of substructure $b$ . Similar, the CCS orientation matrix can be written as + +$$ +\begin{array}{r}{\pmb{\mathsf{E}}_{c,\beta}=\pmb{\mathsf{R}}_{0,b}\left(t,\mathbf{q}\right)\,\pmb{\mathsf{E}}_{c,1,\beta,b}\left(\mathbf{q}_{b}\right)}\end{array} +$$ + +where $\pmb{\mathsf{E}}_{c,1,\beta,b}\;=\;[\pmb{\mathsf{e}}_{1,c,1,\beta,b}\;\pmb{\mathsf{e}}_{2,c,1,\beta,b}\;\pmb{\mathsf{e}}_{3,c,1,\beta,b}]$ is the CCS matrix in the frame of substructure $b$ . Using (2.13) and (2.14) and isolating integrals defined in the substructure frame (marked in blue similar to the substructure components in Section 1.1) using the matrix dot product (E.7), the generalized force (2.11) on and from the aerodynamic force distribution over the substructure $b$ can be written as + +$$ +\begin{array}{l}{{Q_{i,\beta,b}=\mathbf{r}_{0,b,q_{i}}\mathbf{R}_{0,b}\int_{0}^{L_{b}}{\mathbf{f}_{1,\beta,b}\,d s}+\left(\mathbf{R}_{0,b,q_{i}}^{T}\mathbf{R}_{0,b}\right):\int_{0}^{L_{b}}{\left(\mathbf{f}_{1,\beta,b}\mathbf{r}_{a c,1,\beta,b}^{T}+\mathbf{e}_{2,c,1,\beta,b}\mathbf{e}_{1,c,1,\beta,b}^{T}M_{\beta}\right)d s}}\\ {{\ \ \ \ \ \ \ \ \ +\int_{0}^{L_{b}}{\left(\mathbf{r}_{a c,1,\beta,b,q_{i}}^{T}\mathbf{f}_{1,\beta,b}+\mathbf{e}_{1,c,1,\beta,b,q_{i}}^{T}\mathbf{e}_{2,c,1,\beta,b}\,M_{\beta}\right)d s}}\end{array} +$$ + +where $L_{b}$ is the curve­length of blade part on substructure $b$ , the notation $\partial/\partial q_{i}\ =\ ()_{,q_{i}}$ denotes the DOF derivatives, and we have introduced the cross­sectional force vector described in the substructure frame as + +$$ +\mathbf{f}_{1,\beta,b}=\mathbf{E}_{c,1,\beta,b}\mathbf{f}_{c,\beta} +$$ + +Similar, we also introduce a “force vector” in the substructure frame for the subsequent derivations as + +$$ +\begin{array}{r}{\mathbf{m}_{1,\beta,b}=\mathbf{e}_{2,c,1,\beta,b}M_{\beta}}\end{array} +$$ + +whereby the virtual work of the aerodynamic moment on a DOF $q_{i}$ on the substructure $\mathbf{e}_{1,c,1,\beta,b,q_{i}}^{T}\mathbf{e}_{2,c,1,\beta,b}\,M_{\beta}$ can be written as $\mathbf{e}_{1,c,1,\beta,b,q_{i}}^{T}\mathbf{m}_{1,\beta,b}$ , similar to the virtual work of a force. + +The first line of (2.15) are the generalized forces on DOFs translating or rotating the base of the substructure $b$ , and the second line are generalized forces on the DOFs of the substructure itself. Using the superscript $^{\ast}0^{\ast}$ for forces on DOFs moving the substructure base and “1” for forces on substructure DOFs, we split the generalized aerodynamic force into these two categories and rewrite them in the following form: + +$$ +\begin{array}{r l}&{Q_{i,\beta,b}^{0}=\pmb{\Gamma}_{0,b,q_{i}}\pmb{\mathrm{R}}_{0,b}\pmb{\mathrm{f}}_{\beta,b}+\left(\pmb{\mathrm{R}}_{0,b,q_{i}}^{T}\pmb{\mathrm{R}}_{0,b}\right):\pmb{\mathrm{M}}_{\beta,b}}\\ &{Q_{i,\beta,b}^{1}=\displaystyle\int_{0}^{L_{b}}\left(\pmb{\Gamma}_{a c,1,\beta,b,q_{i}}^{T}\pmb{\mathrm{f}}_{1,\beta,b}+\pmb{\mathrm{e}}_{1,c,1,\beta,b,q_{i}}^{T}\pmb{\mathrm{m}}_{1,\beta,b}\right)d s}\end{array} +$$ + +where we have introduced the total aerodynamic forces on the blade ${\mathfrak{f}}_{\beta,b}$ and a ${\mathfrak{3}}{\times}{\mathfrak{3}}$ matrix containing the total aerodynamic moments ${\tt M}_{\beta,b}$ defined as + +$$ +\mathbf{f}_{\beta,b}=\int_{0}^{L_{b}}\mathbf{f}_{1,\beta,b}\,d s\quad\mathrm{and}\quad\mathbf{M}_{\beta,b}=\int_{0}^{L_{b}}\big(\mathbf{f}_{1,\beta,b}\mathbf{r}_{a c,1,\beta,b}^{T}+\mathbf{m}_{1,\beta,b}\mathbf{e}_{1,c,1,\beta,b}^{T}\big)\,d s +$$ + +To show the physical meaning of this total aerodynamic moment matrix, we note that $\mathsf{R}_{0,b,q_{i}}^{T}\bar{\mathsf{R}}_{0,b}$ in (2.18a) is a skew­symmetric matrix with the generic form: + +$$ +\frac{\partial\mathsf{R}^{T}}{\partial q_{i}}\mathsf{R}=\left[\begin{array}{c c c}{0}&{\delta\theta_{z}}&{-\delta\theta_{y}}\\ {-\delta\theta_{z}}&{0}&{\delta\theta_{x}}\\ {\delta\theta_{y}}&{-\delta\theta_{x}}&{0}\end{array}\right] +$$ + +describing the variational rotations $\delta\theta_{x},\,\delta\theta_{y}$ , and $\delta\theta_{x}$ of the substructure base about the three axes in the off­ diagonal elements (follows from Equations (E.1) and (E.5)). If we let $\pmb{r}=\{x,y,z\}^{T}$ and $\pmb{\mathsf{f}}=\{f_{x},f_{y},f_{z}\}^{T}$ be some position and force vectors and derive the moment matrix about the origin of the position vector as + +$$ +{\boldsymbol{\mathsf{M}}}=\mathbf{f}\mathbf{r}^{T}\left[{\begin{array}{l l l}{x f_{x}}&{y f_{x}}&{z f_{x}}\\ {x f_{y}}&{y f_{y}}&{z f_{y}}\\ {x f_{z}}&{y f_{z}}&{z f_{z}}\end{array}}\right] +$$ + +The matrix dot product part of the generalized force on a rotational DOF $q_{i}$ has therefore the generic form + +$$ +\left(\frac{\partial\mathsf{R}^{T}}{\partial q_{i}}\mathsf{R}\right):\mathsf{M}=\delta\theta_{x}\left(y f_{z}-z f_{y}\right)+\delta\theta_{y}\left(z f_{x}-x f_{z}\right)+\delta\theta_{z}\left(x f_{y}-y f_{x}\right) +$$ + +which shows the physical meaning of the second term in (2.18a) as virtual work of the six moment components of the moment matrix (note that its diagonal are not representing moments and therefore obsolete). + +# 2.2 Aerodynamic forces and moment at a cross­section + +This section contains the derivation of the velocity triangles used to compute the aerodynamic forces and moment at an airfoil cross­section along the blade. The local aerodynamic forces defined in the CCS of the cross­section can be computed from the lift and drag forces as + +$$ +f_{x,\beta}=L_{\beta}\sin\alpha_{\beta}-D_{\beta}\cos\alpha_{\beta}\quad{\mathsf{a n d}}\quad f_{y,\beta}=L_{\beta}\cos\alpha_{\beta}+D_{\beta}\sin\alpha_{\beta} +$$ + +where $L_{\beta}$ and $D_{\beta}$ are the lift and drag forces, and $\alpha_{\beta}=$ arctan $(v_{y}/(-v_{x}))$ is the geometrical angle of attack given by the 2D relative flow vector $\{v_{x},v_{y}\}^{T}$ in the $(x,y)$ ­plane of the CCS along blade $\beta$ . The aerodynamic forces and moment about the aerodynamic center are defined as + +$$ +\begin{array}{c}{{\displaystyle{L_{\beta}=\frac{1}{2}\rho c_{\beta}U_{\beta}^{2}C_{L,\beta}^{d y n}}}}\\ {{\displaystyle{D_{\beta}=\frac{1}{2}\rho c_{\beta}U_{\beta}^{2}C_{D,\beta}^{d y n}}}}\\ {{\displaystyle{M_{\beta}=\frac{1}{2}\rho c_{\beta}^{2}U_{\beta}^{2}C_{M,\beta}^{d y n}}}}\end{array} +$$ + +owfh tehree $\rho$ i irfs otilh,e aanird dtehnes iutyn, $c_{\beta}$ aisd tyh ae elroocdayl ncahomridc lceonegftfhic, $U_{\beta}=\sqrt{v_{x}^{2}+v_{y}^{2}}$ eids ,a tiavned f poefe tdh ien tuhnes tpelaandey a ste ients are denot CLd,yβn , $C_{L,\beta}^{d y n},\,C_{D,\beta}^{d y n}$ + +aerodynamics model (cf. Section 2.2.2). Using these definitions of lift and drag and geometric angle of attack, the aerodynamic forces in the $(x,y)$ ­plane of the CCS can be rewritten as + +$$ +f_{x,\beta}=\frac{1}{2}\rho c_{\beta}U_{\beta}\left(v_{y}\,C_{L,\beta}^{d y n}+v_{x}\,C_{D,\beta}^{d y n}\right)\quad\mathrm{~and~}\quad f_{y,\beta}=\frac{1}{2}\rho c_{\beta}U_{\beta}\left(-v_{x}\,C_{L,\beta}^{d y n}+v_{y}\,C_{D,\beta}^{d y n}\right) +$$ + +where the projected 2D relative flow vector $\left\{v_{x},v_{y}\right\}^{T}$ is derived in the next section on the velocity triangles. Note that the lift and drag forces are perpendicular and parallel to the relative flow vector and do not rotate with the airfoil; a rotation of the cross­section due to torsion will only change magnitude of the forces, which only indirectly may change their direction through the dynamic inflow effect on the induced velocities. + +# 2.2.1 Velocity triangles + +The wind field and induced velocity vectors in the ground­fixed inertia frame at a point $p$ on the blade are denoted $\mathbf{w}$ and $\mathbf{v}_{i}$ , and the total relative velocity vector in this frame can be written as + +$$ +\mathbf{v}_{p,r e l,\beta}=\mathbf{w}+\mathbf{v}_{i}-\mathbf{v}_{p,\beta} +$$ + +where the wind field vector $\pmb{\mathsf{w}}=\pmb{\mathsf{w}}\left(t,\pmb{\mathsf{r}}_{p,\beta}\right)$ depends on time and the position of the point $\mathbf{r}_{p,\beta}=\mathbf{r}_{p,\beta}\left(t,\mathbf{q};x,y,s\right)$ and the induced velocity vector $\mathbf{v}_{i}=\mathbf{v}_{i}\left(s,t,\mathbf{q},\dot{\mathbf{q}},\mathbf{x}_{a},\mathbf{w}\right)$ depends on the spanwise position, time (in case of a prescribed rotor speed), and all model state variables. The velocity vector $\mathbf{v}_{p,\beta}=d\mathbf{r}_{p,\beta}/d t$ is the velocity of the point in the inertia frame due to the structural motion of the blade, which we will derive below. The wind field models are described in Appendix D and the wake induction model is described in Section 2.4. + +The relative velocity vector in the CCS can be derived as + +$$ +\begin{array}{r}{\mathbf{\check{v}}_{c,p,r e l,\beta}=\mathbf{\check{E}}_{c,\beta}^{T}\left(\mathbf{\check{w}}+\mathbf{\check{v}}_{i}-\mathbf{\check{v}}_{p,\beta}\right)}\end{array} +$$ + +where $\mathsf{E}_{c,\beta}$ is the orientation matrix of the CCS in the ground­fixed frame. The $x^{.}$ ­ and $y$ ­components of this relative velocity vector are given by + +$$ +\begin{array}{r}{v_{x,c,p,r e l,\beta}=\pmb{\mathbf{e}}_{1,c,\beta}^{T}\left(\pmb{\mathbf{w}}+\pmb{\mathbf{v}}_{i}-\pmb{\mathbf{v}}_{p,\beta}\right)\quad\mathrm{~and~}\quad v_{y,c,p,r e l,\beta}=\pmb{\mathbf{e}}_{2,c,\beta}^{T}\left(\pmb{\mathbf{w}}+\pmb{\mathbf{v}}_{i}-\pmb{\mathbf{v}}_{p,\beta}\right)}\end{array} +$$ + +where $\mathfrak{e}_{1,c,\beta}$ and $\mathbf{e}_{2,c,\beta}$ are the $x^{.}$ ­ and $y$ ­unit­vectors of the CCS in the ground­fixed frame. + +Figure 2.2 illustrates the two velocity triangles at two different inplane locations of the CCS needed for the computation of the aerodynamic forces and moment in CASEStab. The velocity triangle at the torsional point (TP) is used to compute the magnitude $U_{\beta}$ and direction relative to the chord $\alpha_{\beta}$ of the relative flow vector. The TP is the point where the translations of the cross­section is defined in the kinematic model of the substructure (the origin of the reference frame in case of the co­rotational beam model), i.e., a rotation of the cross­section will not lead to any translations of this point. In contrast to $\alpha_{\beta}$ , the geometric angle of attack $\alpha_{3/4,\beta}$ used to compute the circulatory lift must be based on the velocity triangle set up for the collocation point at the three­quarter chord point [7], which will include downwash due to rate of torsional blade deformation. + +The relative flow speed and the two geometric angles of attack are computed as + +$$ +\begin{array}{r l}&{\iota=\sqrt{v_{x,c,T P,r e t,\beta}^{2}+v_{y,c,T P,r e t,\beta}^{2}}\ ,\quad\alpha_{\beta}=\mathsf{a r c t a n}\,\frac{v_{y,c,T P,r e t,\beta}}{-v_{x,c,T P,r e t,\beta}}\ \ \mathsf{a n d}\ \ \alpha_{3/4,\beta}=\mathsf{a r c t a n}\,\frac{v_{y,c,C P,r e t,\beta}}{-v_{x,c,C P,r e t,\beta}}\ }\end{array} +$$ + +where the relative flow velocities are defined at $p=T P$ and $p=C P$ using (2.25). + +# Relative flow vector on a substructure of the blade + +The blade deformation may be described by the individual deformation and base motion of several substruc­ tures. We write the position of a blade point $p$ in a form similar to the structural kinematics Eq. (1.12): + +$$ +\begin{array}{r}{\mathbf{r}_{p,\beta,b}\left(t,\mathbf{q};x,y,s\right)=\mathbf{r}_{0,b}\left(t,\mathbf{q}\right)+\mathbf{R}_{0,b}\left(t,\mathbf{q}\right)\mathbf{r}_{p,1,\beta,b}\left(\mathbf{q}_{b};x,y,s\right)}\end{array} +$$ + +where the translations $\boldsymbol{\mathsf{r}}_{0,b}$ and rotations $\pmb{\mathsf{R}}_{0,b}$ of the substructure number $b$ are given by (1.13) and they may depend explicit on time $t$ in case of prescribed rotor speeds. The local position vector in the substructure + +![](images/9f4bef8a8ad3880f61f4b8e644a9e5035838d6dd9f9cf870563a547bc91c921a.jpg) +Figure 2.2: Velocity triangle at the torsional point (TP) used for computing the relative flow speed $U_{\beta}$ and the geometric angle of attack $\alpha_{\beta}$ that defines the orientations of the lift $L_{\beta}$ and drag $D_{\beta}$ forces acting at the aerodynamic center (AC). Another velocity triangle is defined at the collocation point (CP) to compute the geometric angle of attack $\alpha_{3/4,\beta}$ including downwash from torsional rate in the circulatory lift modeling. + +frame $\boldsymbol{\mathsf{r}}_{p,1,\beta,b}$ is a function of the DOFs on the blade $\mathbf{q}_{b}$ , the cross­sectional coordinates $x,y$ , and the spanwise coordinate $s$ on the substructure. We assume, as in the structural kinematics, that the deflection of the blade has no explicit dependency of time. The velocity vector can be derived in the ground­fixed inertia frame as + +$$ +\mathbf{\mu}_{p,\beta,b}=\frac{d\mathbf{r}_{p,\beta,b}}{d t}=\dot{\mathbf{r}}_{0,b}+\sum_{i=1}^{N_{D}}\mathbf{r}_{0,b,q_{i}}\dot{q}_{i}+\dot{\mathbf{R}}_{0,b}\mathbf{r}_{p,1,\beta,b}+\left(\sum_{i=1}^{N_{D}}\mathbf{R}_{0,b,q_{i}}\dot{q}_{i}\right)\mathbf{r}_{p,1,\beta,b}+\mathbf{R}_{0,b}\left(\sum_{i=1}^{N_{D}}\mathbf{r}_{p,1,\beta,b,q_{i}}\dot{q}_{i}\right) +$$ + +and we write the CCS orientation matrix in ground­fixed frame as + +$$ +\begin{array}{r}{\pmb{\mathsf{E}}_{c,\beta}=\pmb{\mathsf{R}}_{0,b}\left(t,\mathbf{q}\right)\,\pmb{\mathsf{E}}_{c,1,\beta,b}\left(\mathbf{q}_{b}\right)}\end{array} +$$ + +where $\pmb{\mathsf{E}}_{c,1,\beta,b}\:=\:[\pmb{\mathsf{e}}_{1,c,1,\beta,b}\:\pmb{\mathsf{e}}_{2,c,1,\beta,b}\:\pmb{\mathsf{e}}_{3,c,1,\beta,b}]$ is the CCS matrix in the substructure frame. Thus, the $x^{.}$ ­ and $y$ ­components of the relative velocity vector in the CCS can be computed as + +
T RO,b T RT,b RT,b ro,b Ux,c,p,rel,β,b w+ND RO,b ro,b,qi Qi 2=1 ND(2.30a)
ND RT,b Ro,b T RO,b p,1,β,be 1,c,1,β,bRo,b,qi Qi T
=1 RT Vy,c,p,rel,β,b e w+ T e.L ro,bT M i=1 NDIp,1,β,b,qi Qi
M RO,b bro,b,qiQi(2.30b)
2,c,1,β,b 2,c,1,β,b ND RT,b Ro,b (rp,1,β,be2,c,1,β,b) T i=10,b (p,1,β,be2,c,1,β,b, 1i=1 ND e2,c,1,β,b T M p,1,β,b,qiQi i=1
+ +where we have blue marked the vectors that are defined in the substructure frame and used the matrix dot product (E.7) to isolate these components of the local 2D flow vector. Note that only the last term in these local velocity components describe the contribution from the deformation velocities of the substructure $b$ on blade $\beta$ . + +# 2.2.2 Unsteady aerodynamic model + +The dynamic model of the unsteady airfoil aerodynamics is very similar to the Beddoes­Leishman type model described in [5], except that the attached flow lift curve is not presumed to be linear but follows the user­defined + +lift curve until the angles of attack where flow separation initiates, and except that the part of the unsteadiness on the moment coefficient is omitted. + +# Stationary steady­state airfoil aerodynamics + +For computation of the aerodynamic power and the stationary steady­state deflection of the blades in a uniform wind field and no gravity or tower shadow, the aerodynamic coefficients can be obtained from the user­defined airfoil polar (denoted by superscript $^{s t}$ ) as + +$$ +\begin{array}{r l}&{C_{L,\beta}^{d y n}=C_{L,\beta}^{s t}\left(\alpha_{\beta}\right)}\\ &{C_{D,\beta}^{d y n}=C_{D,\beta}^{s t}\left(\alpha_{\beta}\right)}\\ &{C_{M,\beta}^{d y n}=C_{M,\beta}^{s t}\left(\alpha_{\beta}\right)}\end{array} +$$ + +where $\alpha_{\beta}$ is the geometric angle of attack given by Equation (2.26). Note that in the stationary steady­state the two geometric angles of attack in (2.26) is identical. + +# Unsteady airfoil aerodynamics + +The unsteady aerodynamic model can be written as [5] + +# 2.3 Aerodynamic discretization of a blade + +In CASEStab, the aerodynamic forces and moment are computed at $N_{a}$ aerodynamic calculation points (ACPs) along the blade curve length and linear interpolation is applied between these sections. This aerodynamic discretization is defined by the user independent of the structural discretization to ensure that the variations of aerodynamic forces in the root sections and at the blade tip, as well as variations caused by aerodynamic devices (e.g. vortex generators) are captured with sufficient accuracy. + +The linear interpolation of forces and moment and the integrations over the blade are defined for the blade curve­coordinate $s$ ; however, this curve­coordinate is not well­defined for the ACPs because they will depend on the blade data resolution defined by the user. Eventually, the integrations over the blade must be sub­divided into integrations over each substructure $b=b_{\beta},\ldots,b_{\beta}+B_{\beta}-1$ of the blade, and possibly further sub­divided into integrations over each element of each substructure (e.g. for the co­rotational beam formulation). As this structural resolution is independent of the aerodynamic resolution, the definition of a curve­coordinate could be inconsistent. Thus, to link the spanwise positions of the $N_{a}$ ACPs on the blade to positions on the substructure, we define them by the coordinates $z_{a,j}\;(j=1,2,.\,.\,.\,,N_{a})$ along the axis perpendicular to the flange of the blade, which is the same of the pitch axis for a conventional blade mounted on a pitch bearing. Note that $z_{a,1}=0$ and the coordinate of the blade tip is smaller than the blade curve­length $z_{a,N_{a}}2.5}\end{array}\right. +$$ + +# (2.32) + +# 2.4.1 Stationary steady­state rotor aerodynamics + +In case of uniform wind field (without shears, yaw, upflow, or turbulence), exclusion of tower shadow effects, and no periodic forcing of the blades (no gravity or control actions), the rotor + +# 2.4.2 Periodic steady­state rotor aerodynamics + +# 3 Aeroelastic equations of motion + +The structural motion and the aerodynamic forces described in the previous two chapters are now coupled in a closed set of aeroelastic equations of motion. The generalized aerodynamic forces on the structural DOFs are given on a generic form in Equation (2.18) for the integration of the aerodynamic forces and moment distribution over a single substructure. The relative flow at the airfoil section of a blade depend on the structural motion as described on generic form in Section 2.2. The aerodynamic forces and moment at the airfoil section depend on this relative flow, the wind field, the induced velocities, and local dynamic effects of shed vorticity and dynamic stall. The aerodynamic discretization of the blade into aerodynamic calculation points (ACPs) and the piecewise linear functions used to describe the distributions of forces and moment are introduced in Section 2.3. The unsteady airfoil aerodynamic model is described in Section 2.2.2. The wake model of dynamic inflow (the unsteadiness of the induced velocities) depend on the local thrust and torque coefficients and the yaw and upflow angles as described in Section 2.4. + +The chapter is structured as follows: Section 3.1 + +# 3.1 Generalized aerodynamic forces for different substructure types + +The total generalized aerodynamic force from a blade on each structural DOF is given in (2.12) by the sum of the contributions from aerodynamic forces and moment distributions on each substructure $b=b_{\beta},\ldots,b_{\beta}+B_{\beta}-1$ of the $B_{\beta}$ substructures comprising blade number $\beta$ . These contributions are given by the integrals over each substructure in Equations (2.18) and (2.19). In this section, we derive expressions for these integrals for the different types of substructures in CASEStab using piecewise linear functions of the aerodynamic forces and moment distributions over the substructure. The sectional force and moment vectors at the ACPs are denoted $\mathfrak{f}_{1,\beta,b,j}$ and $\mathbf{m}_{1,\beta,b,j}$ , where the first subscript “1” refers to the fact that the vector has been computed in or transformed to the frame of the substructure number $b$ being the third subscript, the second subscript $\beta$ is the blade number, and the fourth and last subscript is the ACP number $j=1,\dots,N_{a}$ on the blade. + +We define an index vector $j_{b,k}$ for $k=1,\ldots,N_{a,b}$ containing the indices of the $N_{a,b}$ ACPs which sectional forces and moments define the forces and moment distributions over the substructure $b$ . We collect these forces and moments needed for the linear interpolation over substructure number $b$ in the following vectors: + +$$ +\mathbf{f}_{a l l,1,\beta,b}=\left\{\begin{array}{c}{\mathbf{R}_{0,b}^{T}\mathbf{R}_{0,b-1}\mathbf{f}_{1,\beta,b-1,j_{b,1}}}\\ {\mathbf{f}_{1,\beta,b,j_{b,2}}}\\ {\mathbf{f}_{1,\beta,b,j_{b,3}}}\\ {\vdots}\\ {\mathbf{f}_{1,\beta,b,j_{b,N_{a,b}-1}}}\\ {\mathbf{R}_{0,b}^{T}\mathbf{R}_{0,b+1}\mathbf{f}_{1,\beta,b+1,j_{b,N_{a,b}}}}\end{array}\right\}\;\;\mathrm{and}\;\;\mathbf{m}_{a l l,1,\beta,b}=\left\{\begin{array}{c}{\mathbf{R}_{0,b}^{T}\mathbf{R}_{0,b-1}\mathbf{m}_{1,\beta,b-1,j_{b,1}}}\\ {\mathbf{m}_{1,\beta,j_{b,2}}}\\ {\mathbf{m}_{1,\beta,b,j_{b,3}}}\\ {\vdots}\\ {\mathbf{m}_{1,\beta,b,j_{b,N_{a,b}-1}}}\\ {\mathbf{R}_{0,b}^{T}\mathbf{R}_{0,b+1}\mathbf{m}_{1,\beta,b+1,j_{b,N_{a,b}}}}\end{array}\right\} +$$ + +where the first and last vectors defined on adjacent substructures must be transformed into the frame of sub­ structure $b$ by the orientation matrices $\mathbf{R}_{0,b-1}$ , $\pmb{\mathsf{R}}_{0,b}$ , and $\mathsf{R}_{0,b+1}$ for each involved substructure (1.13). If the blade is described by a single substructure then $j_{b}=1$ and $N_{a,b}=N_{a}$ , and the above vectors reduce to + +$$ +\mathbf{f}_{a l l,1,\beta,b}=\left\{\begin{array}{l}{\mathbf{f}_{1,\beta,b,j_{b,1}}}\\ {\mathbf{f}_{1,\beta,b,j_{b,2}}}\\ {\qquad\vdots}\\ {\mathbf{f}_{1,\beta,b,j_{b,N_{a}}}}\end{array}\right\}\,\,\,\mathbf{and}\,\,\,\mathbf{m}_{a l l,1,\beta,b}=\left\{\begin{array}{l}{\mathbf{m}_{1,\beta,b,j_{b,1}}}\\ {\mathbf{m}_{1,\beta,b,j_{b,2}}}\\ {\qquad\vdots}\\ {\mathbf{m}_{1,\beta,b,j_{b,N_{a}}}}\end{array}\right\} +$$ + +as one special case of (3.1). The other two special cases are when $b$ is the first or last of several substructures on the blade. + +The outcome of this section are matrices that transform the $N_{a}$ forces and $N_{a}$ moments of the ACPs on the blade to the total aerodynamic force and moment, as well as the + +Note that we omit the subscript $\beta$ for the blade number in the following sections because we are referring to a single blade. + +# 3.1.1 Rigid body substructure + +To be derived ... + +# 3.1.2 Co­rotational beam substructure + +Figure 3.1 shows a 2D illustration of the aerodynamic force distribution over a blade and an element of a co­rotational beam element substructure of the blade. The blade consists of two substructures $b-1$ and $b$ described by several beam elements with their end nodes are marked by $\bullet$ . The force distribution over the blade is a piecewise linear function with the $N_{a}$ ACPs (marked by ◦) as break­points. The piecewise linear function over the element $n$ of substructure $b$ is here defined by four ACPs, written as $N_{a,b,n}\,=\,4$ , i.e., the force (and moment) distribution over this element is described by two internal break­points and two end­break­ points on the adjacent elements. We define a vector $\zeta_{b,n,k}$ for $k=1,\ldots,N_{a,b,n}$ containing the non­dimensional element coordinate of these break­points for each element $n$ . Note that the end­break­points will always lie at the element nodes $\left\langle\zeta\right.=\pm1)$ ) or outside the element, i.e., $\zeta_{b,n,1}\,\leq\,-1$ and $\zeta_{b,n,N_{a,b,n}}\,\geq\,1$ . In case that these points lie on the adjacent elements, their element coordinates can be computed as + +$$ +\zeta_{b,n,1}=-1-\frac{l_{n-1}}{l_{n}}\left(1-\zeta_{b,n-1,N_{a,b,n-1}-1}\right)\quad\mathrm{~and~}\quad\zeta_{b,n,N_{a,b,n}}=1+\frac{l_{n+1}}{l_{n}}\left(1+\zeta_{b,n+1,2}\right) +$$ + +where $l_{n}$ is the initial lengths of elements $n\,=\,1,\ldots,N_{e,b}$ . Note that special rules apply to the first and last elements of a substructure, where the adjacent element may lie on an adjacent substructure. + +A curve integral over the co­rotational beam substructure $b$ is the sum of integrals over each element $n$ as + +$$ +\int_{0}^{L_{b}}\left(\mathbf{\Phi}\right)d\sigma=\sum_{n=1}^{N_{e,b}}\frac{l_{n}}{2}\int_{-1}^{1}\left(\mathbf{\Phi}\right)d\zeta=\sum_{n=1}^{N_{e,b}}\frac{l_{n}}{2}\left(\sum_{m=1}^{N_{a,b,n}-1}\int_{a_{b,n,m}}^{b_{b,n,m}}\left(\mathbf{\Phi}\right)d\zeta\right) +$$ + +where each curve integral over the element length is a line integral over the non­dimensional element coordinate $\zeta$ between element its nodes at $\zeta\,=\,\pm1$ . To ensure that the integral kernels will be single polynomials (not piecewise), the line integral over each element is further sub­divided into line integrals between the nodes and the $N_{a,b,n}-2$ internal break­points on the element. The limits of these sub­divided integrals are + +$$ +\begin{array}{r c l}{{(a_{b,n,1},b_{b,n,1})}}&{{=}}&{{(-1,\zeta_{b,n,2})}}\\ {{(a_{b,n,2},b_{b,n,2})}}&{{=}}&{{(\zeta_{b,n,2},\zeta_{b,n,3})}}\\ {{}}&{{\vdots}}&{{}}\\ {{(a_{b,n,N_{a,b,n}-1},b_{b,n,N_{a,b,n}-1})}}&{{=}}&{{(\zeta_{b,N_{a,b,n}-1},1)}}\end{array} +$$ + +Note that in case there are no ACPs on an element, when $N_{a,b,n}=2$ , only a single integration over the element is needed $\left(a_{b,n,1},b_{b,n,1}\right)\,=\,\left(-1,1\right)$ . The piecewise linear distributions of the aerodynamic force and moment vectors over the element can be written on generic forms as + +$$ +\mathbf{f}_{1,b,n}\left(\zeta\right)=\left\{\begin{array}{c c}{\sum_{r=0}^{1}\mathbf{f}_{1,b,n,0,r}\zeta^{r}}&{a_{b,n,1}\leq\zeta
1 z [m]Distance tosection frombase ofsubstructure(bladeroot flange)along the third (pitch)
2 xref [m]axis of the substructure (blade) frame. Primary (edge) axis distance from third substructure (pitch) axis to the origin (the
3 Yref[m]Secondary (flap) axis distance from third substructure (pitch) axis to the origin (the
4 Oref[deg]reference point) of the cross-sectional reference frame. Eachelementcoordinatesystemisrotated theelement-averageofthisanglefrom an initial orientation where its secondary (flap) axis has a zero primary (edge) axis component in the substructure (blade) frame.
+ +# Mass properties: + +
5 m Xcg[kg/m] Massperunit-lengthofthesubstructure.Distance on primary axis of the reference frame from reference point to mass center
6[m](CG).
7ycg [m]Distance on secondary axis of the reference frame from reference pointto mass center(CG).
8rc [m]Radius gyration for rotation about the primary inertia axis with CG as origin.
9ry [m]Radius gyration for rotation about the secondary inertia axis with CG as origin.
10OrC [deg]Angular offset of the primary inertia axis from the primary axis of the reference frame.
+ +# Isotropic section: + +Stiffness properties (isotropic section $=13$ columns or full 6x6 compliance matrix $=$ 21 columns) + + +
11 αea [m]
12 yea [m]of bending (elastic axis).
13 sc [m]of bending (elastic axis).
14 ysc [m]center. Distance on secondary axis of the reference frame from reference point to the shear
15 Obend [deg]center.
16 E [Pa](structural twist). Average elastic Young's modulus of the cross-section.
17 G [Pa]Average shear modulus of the cross-section.
18 A [m²]Areaofcross-section.
19 Ix [m4]Moment of inertia for bending about the primary bending axis.
20 Iy [m4]Moment of inertia for bending about the secondary bending axis.
21 K [m4]Momentofinertiafortorsion.
22 kx[-]Shear correctionfactor in the direction of theprimaryreference axis.
23 ky [-]Shear correction factor in the direction of the secondary reference axis.
+ +# Full 6x6 compliance matrix: + +1­31 C [SI] Upper triangle 21 elements of the symmetric $6\!\times\!6$ compliance matrix row by row $[C_{11},C_{12},.\,.\,.\,,C_{16},C_{22},.\,.\,.\,,C_{66}]$ given in the reference frame. + +# B Aerodynamic data input and process­ ing + +# B.1 Input files + +
1 ≈[m]Distance to section from base of substructure (blade root flange) along the third (pitch)
2Cref [m]axisof thesubstructure(blade)frame. Primary (edge) axis distance from third substructure (pitch) axis to the reference point
3Yref [m]of thecross-sectional referenceframe. Secondary (flap) axis distance from third substructure (pitch) axis to the reference
4-6Px,Py,Φz [deg]pointofthecross-sectionalreferenceframe. Pseudo vector of finite rotation of Chord Coordinate System (CCS) from the substruc- ture (blade) frame. If Φx = Φy = O, then the normal to the airfoil plane is defined by the tangent of the reference point curve and the inplane axes of the CCS are defined
7such that its secondary axis has zero primary axis component in the blade frame before rotating them about the tangent by Φz.
8 9c[m] tret [%]Chordlengthofairfoilordiameterofcircle. Relative thickness of airfoil (100% = circle).
ac [m]Distance on primary axis of the reference frame from reference point to the aerody-
10yac [m]namic center (AC). Distance on secondary axis of the reference frame from reference point to the aero-
11αac [-]dynamic center (AC).
12ipc[-]Numberof theairfoilsetinHAWC2typepolardatafiles.
+ +# C Dynamic stall model data + +This appendix contains description of the pre­processing of the data for the dynamic stall model based on the user­defined airfoil polars. + +# D Wind field + +# E Mathematics + +This appendix contains different mathematical formulas and derivations that are used in the theory manual. + +# E.1 Definitions + +Let ${\mathfrak{s}}\{\}$ be a skew­symmetric operator converting a vector to a skew­symmetric matrix defined as + +$$ +\begin{array}{r}{\pmb{\ S}\{\omega\}=\left[\begin{array}{c c c}{0}&{-\omega_{3}}&{\omega_{2}}\\ {\omega_{3}}&{0}&{-\omega_{1}}\\ {-\omega_{2}}&{-\omega_{1}}&{0}\end{array}\right]}\end{array} +$$ + +where $\pmb{\omega}=\{\omega_{1},\omega_{2},\omega_{3}\}^{T}$ is any vector. If S is a skew­symmetric matrix and $\blacktriangle$ is a arbitrary square matrix of the same size then $\mathbb{A}^{T}\mathbb{S}\mathbb{A}$ is also a skew­symmetric matrix. + +# E.2 Finite Rotations + +Any finite rotation can be described by a unit­vector $\pmb{\mathsf{n}}=\{n_{x},n_{y},n_{z}\}^{T}$ and an angle $\phi$ the rotation matrix: + +$$ +\boldsymbol{\mathsf{R}}=\boldsymbol{\mathsf{I}}+\boldsymbol{\mathsf{N}}\sin\phi+\boldsymbol{\mathsf{N}}^{2}\left(1-\cos\phi\right) +$$ + +where + +$$ +\mathsf{\pmb{N}}=\left[\begin{array}{c c c}{0}&{-n_{z}}&{n_{y}}\\ {n_{z}}&{0}&{-n_{x}}\\ {-n_{y}}&{n_{x}}&{0}\end{array}\right]=\mathsf{\pmb{S}}\{\mathsf{n}\} +$$ + +is a skew­symmetric matrix. Using the unity of the rotation vector length $\|\pmb{\mathsf{n}}\|=1$ , it can be shown that + +$$ +\mathsf{N}^{3}=-\mathsf{N} +$$ + +and using this property, the derivative of the rotation matrix with respect to the angle for a fixed rotation vector becomes + +$$ +\frac{\partial\mathsf{R}}{\partial\phi}=\mathsf{N}\cos\phi+\mathsf{N}^{2}\sin\phi=\mathsf{R}\mathsf{N} +$$ + +and therefore + +$$ +\frac{\partial^{2}\mathsf{R}}{\partial\phi^{2}}=\mathsf{R N}^{2} +$$ + +# E.3 Miscellaneous formulas + +$$ +\mathbf{u}^{T}\mathbf{U}^{T}\mathbf{V}\mathbf{v}=\mathsf{T r}\left(\left(\mathbf{U}^{T}\mathbf{V}\right)\left(\mathbf{v}\mathbf{u}^{T}\right)\right)\equiv\left(\mathbf{U}^{T}\mathbf{V}\right):\left(\mathbf{v}\mathbf{u}^{T}\right) +$$ + +![](images/52979061bffa1fc4421a09d98e07e46c7a71bffe08fdcfb0be650f5a67be09a8.jpg) + +University of Southern Denmark + +Indsæt co-branding logo 2 + +Phone: +45 6550 1000 +sdu@sdu.dk +www.sdu.dk \ No newline at end of file diff --git a/书籍/力学书籍/CASEstab_theory_manual/auto/CASEstab_theory_manual_origin.pdf b/书籍/力学书籍/CASEstab_theory_manual/auto/CASEstab_theory_manual_origin.pdf new file mode 100644 index 0000000..2ab4917 Binary files /dev/null and b/书籍/力学书籍/CASEstab_theory_manual/auto/CASEstab_theory_manual_origin.pdf differ diff --git a/书籍/力学书籍/CASEstab_theory_manual/auto/images/0db388bbd1194bf9846245e1fe9ae46150fbaacce964eab7515591c209fd7103.jpg b/书籍/力学书籍/CASEstab_theory_manual/auto/images/0db388bbd1194bf9846245e1fe9ae46150fbaacce964eab7515591c209fd7103.jpg new file mode 100644 index 0000000..96076bf Binary files /dev/null and b/书籍/力学书籍/CASEstab_theory_manual/auto/images/0db388bbd1194bf9846245e1fe9ae46150fbaacce964eab7515591c209fd7103.jpg differ diff --git a/书籍/力学书籍/CASEstab_theory_manual/auto/images/486f1b76093ec6a8ff87122bde054e63ad3a2e764e44c24f8f7b740bbff29d47.jpg b/书籍/力学书籍/CASEstab_theory_manual/auto/images/486f1b76093ec6a8ff87122bde054e63ad3a2e764e44c24f8f7b740bbff29d47.jpg new file mode 100644 index 0000000..b8e2cae Binary files /dev/null and b/书籍/力学书籍/CASEstab_theory_manual/auto/images/486f1b76093ec6a8ff87122bde054e63ad3a2e764e44c24f8f7b740bbff29d47.jpg differ diff --git a/书籍/力学书籍/CASEstab_theory_manual/auto/images/52979061bffa1fc4421a09d98e07e46c7a71bffe08fdcfb0be650f5a67be09a8.jpg b/书籍/力学书籍/CASEstab_theory_manual/auto/images/52979061bffa1fc4421a09d98e07e46c7a71bffe08fdcfb0be650f5a67be09a8.jpg new file mode 100644 index 0000000..2f5ae8f Binary files /dev/null and b/书籍/力学书籍/CASEstab_theory_manual/auto/images/52979061bffa1fc4421a09d98e07e46c7a71bffe08fdcfb0be650f5a67be09a8.jpg differ diff --git a/书籍/力学书籍/CASEstab_theory_manual/auto/images/662c7eec9fa38e47ccb8876191c81b95627aefb4688ddfcaf96244ee7834a7b7.jpg b/书籍/力学书籍/CASEstab_theory_manual/auto/images/662c7eec9fa38e47ccb8876191c81b95627aefb4688ddfcaf96244ee7834a7b7.jpg new file mode 100644 index 0000000..e459e2e Binary files /dev/null and b/书籍/力学书籍/CASEstab_theory_manual/auto/images/662c7eec9fa38e47ccb8876191c81b95627aefb4688ddfcaf96244ee7834a7b7.jpg differ diff --git a/书籍/力学书籍/CASEstab_theory_manual/auto/images/7b9d537e5e701f6c5058df02575854b6684a5dbbefd131886d1374e7de03a5f6.jpg b/书籍/力学书籍/CASEstab_theory_manual/auto/images/7b9d537e5e701f6c5058df02575854b6684a5dbbefd131886d1374e7de03a5f6.jpg new file mode 100644 index 0000000..e840dfc Binary files /dev/null and b/书籍/力学书籍/CASEstab_theory_manual/auto/images/7b9d537e5e701f6c5058df02575854b6684a5dbbefd131886d1374e7de03a5f6.jpg differ diff --git a/书籍/力学书籍/CASEstab_theory_manual/auto/images/9f4bef8a8ad3880f61f4b8e644a9e5035838d6dd9f9cf870563a547bc91c921a.jpg b/书籍/力学书籍/CASEstab_theory_manual/auto/images/9f4bef8a8ad3880f61f4b8e644a9e5035838d6dd9f9cf870563a547bc91c921a.jpg new file mode 100644 index 0000000..9bfa339 Binary files /dev/null and b/书籍/力学书籍/CASEstab_theory_manual/auto/images/9f4bef8a8ad3880f61f4b8e644a9e5035838d6dd9f9cf870563a547bc91c921a.jpg differ diff --git a/书籍/力学书籍/CASEstab_theory_manual/auto/images/ac7e6ec1853ac26e03566d17b33a6f52f65da377c47d049e52eeb6ca59a149a3.jpg b/书籍/力学书籍/CASEstab_theory_manual/auto/images/ac7e6ec1853ac26e03566d17b33a6f52f65da377c47d049e52eeb6ca59a149a3.jpg new file mode 100644 index 0000000..49d3f10 Binary files /dev/null and b/书籍/力学书籍/CASEstab_theory_manual/auto/images/ac7e6ec1853ac26e03566d17b33a6f52f65da377c47d049e52eeb6ca59a149a3.jpg differ diff --git a/书籍/力学书籍/CASEstab_theory_manual/auto/images/b5ba4db7f3a71c7da780405382854e8e4799f67fdcd07b88d18c6e9f3aeb3d74.jpg b/书籍/力学书籍/CASEstab_theory_manual/auto/images/b5ba4db7f3a71c7da780405382854e8e4799f67fdcd07b88d18c6e9f3aeb3d74.jpg new file mode 100644 index 0000000..c1ba48d Binary files /dev/null and b/书籍/力学书籍/CASEstab_theory_manual/auto/images/b5ba4db7f3a71c7da780405382854e8e4799f67fdcd07b88d18c6e9f3aeb3d74.jpg differ diff --git a/书籍/力学书籍/CASEstab_theory_manual/auto/images/ded3773ba1574335e77bb328c71ac085f65704f51d645beddfb0bcb26f003b39.jpg b/书籍/力学书籍/CASEstab_theory_manual/auto/images/ded3773ba1574335e77bb328c71ac085f65704f51d645beddfb0bcb26f003b39.jpg new file mode 100644 index 0000000..08efaae Binary files /dev/null and b/书籍/力学书籍/CASEstab_theory_manual/auto/images/ded3773ba1574335e77bb328c71ac085f65704f51d645beddfb0bcb26f003b39.jpg differ diff --git a/书籍/力学书籍/CASEstab_theory_manual/auto/images/f0dbd623f7db9b0f64e64930ae8188cb0b48043f7c41d9f97876880f735c4b3f.jpg b/书籍/力学书籍/CASEstab_theory_manual/auto/images/f0dbd623f7db9b0f64e64930ae8188cb0b48043f7c41d9f97876880f735c4b3f.jpg new file mode 100644 index 0000000..af9ee2e Binary files /dev/null and b/书籍/力学书籍/CASEstab_theory_manual/auto/images/f0dbd623f7db9b0f64e64930ae8188cb0b48043f7c41d9f97876880f735c4b3f.jpg differ diff --git a/书籍/力学书籍/input/CASEstab_theory_manual.pdf b/书籍/力学书籍/input/CASEstab_theory_manual.pdf new file mode 100644 index 0000000..b91c9a5 Binary files /dev/null and b/书籍/力学书籍/input/CASEstab_theory_manual.pdf differ