| Module | States | Inputs | Outputs |
| InflowWind (IfW) | None | Positions where the undisturbed (inflow) wind will be indino Disturbances of horizontal wind speed, power-law shear exponent,andwind-propagation direction | Undisturbed (inflow) wind velocity at input positions User-selected wind-inflow outputs |
| ServoDyn (SrvD) | None | Nacelle-yaw angle and rate Generator speed | Blade-pitch-angle command (independent) Nacelle-yaw moment Generator torque and electrical power User-selected control and electrical-drive outputs Translational displacements, orientations, |
| ElastoDyn (ED) | Position and orientation and translational and rotational velocities of the floating platform (continuous states) Relative bending-mode degrees of freedom of the bladesandtowerandtheir first-time derivatives (continuous states) Relative rotational displacement and velocity off the nacelle-yaw, generator azimuth, shaft torsion, and rotor-teeter (continuous states) | Applied point forces and moments distributed along the blades and tower Applied point forces and moments lumped on the hub, nacelle, and platform Blade-pitch-angle command (both independent and rotor- collective) Nacelle-yaw moment Generator torque | translational androtationalvelocities,and translational and rotation accelerations of analysis nodes along the blades and tower Translational displacements, orientations, translational androtationalvelocities,and translationalandrotationaccelerationsoftheblade- root, nacelle, and platform reference points Translational displacement,orientation,and rotational velocity of the hub reference point Nacelle-yaw angle and rate Generator speed User-selected structural outputs |
| BeamDyn (BD) | Positions and orientations and translational and rotational velocitiesofnodes distributed along eachblade (continuous states) | Translational displacements, orientations, translational androtationalvelocities,andtranslationalandrotational accelerations of eachbladeroot Applied point forces and moments distributed along each blade Applied line (per-unit length) forces and moments distributedalongeachblade | Reaction point forces and moments lumped at each bladeroot Translational displacements,orientations, translational and rotational velocities, and translationalandrotationaccelerationsofnodes along each blade User-selected structural outputs |
| AeroDyn (AD) | Inflow angle at each blade analysis node/airfoil (constraint states) | Translational displacements, orientations, and translational velocities of analysis nodes along the blades and tower Orientation of the blade-root reference point Translational displacement, orientation, and rotational velocity of the hub reference point Undisturbed (inflow) wind velocities at analysis nodes along the blades and tower | Aerodynamic-applied line(per-unit length) forces and moments distributed along the blades and tower User-selectedaerodynamicoutputs |
| HydroDyn (HD) | State-space-based wave- excitation states (continuous states) State-space-based wave- radiation states (continuous states) | Translational displacements,orientations, translational and rotational velocities, and translational and rotational accelerations of analysis nodes distributed along the floating platform Translational displacement,orientation,translationand rotational velocities, and translational and rotational accelerations of the platform reference point Disturbanceofwaveelevationattheplatformreference point | Hydrodynamic-applied point forces and moments distributed along the floating platform Hydrodynamic-applied line (per-unit length) forces and moments distributed along the floating platform Hydrodynamic-applied point force and moment at the platform reference point User-selected hydrodynamic outputs |
| MAP++ (MAP) | Horizontal and vertical tensions at thefairleadof each mooring line (constraint states) Positions of each connect node (constraint states for | Translational displacements of each fairlead | Reaction point forces (tensions) lumped at each fairlead User-selected mooring outputs |
+
+The linearization of the HydroDyn (HD) module applies to both the strip-theory solution, potential-flow solution, or a hybrid combination of the two. Linear state-space-based waveexcitation and wave-radiation models have been added to HydroDyn to enable linearization of the potential-flow solution, including wave-excitation loads (with diffraction) from disturbances of wave elevation and wave-radiation loads with free-surface memory effects.
+
+The linear state-space-based wave-radiation model was previously available in HydroDyn within FAST v8. As shown in [5], the model approximates the convolution term in the Cummins equation as shown in Eq. (1), with a linear statespace model as shown in Eq. (2), where $F_{R d t n}$ is the waveradiation loads; $K_{R d t n}$ is the wave-radiation-retardation kernels; $\dot{q}_{P t f m}$ is the floating platform translation and rotational velocities; $x_{R d t n}$ are the wave-radiation states; and $A_{R d t n}$ ,
+
+$B_{R d t n}$ , and $C_{R d t n}$ are the matrices of the linear state-space wave-radiation model. The more wave-radiation states there are, the better the state-space model can approximate the convolution. Within HydroDyn, the user can select either the convolution-based wave-radiation model, which is implemented via numerical convolution using discrete-time states (which, if linearized, would overly complicate application of the linear state-space model), or the linear statespace-based wave-radiation model, which is implemented using matrices derived from the SS_Fitting preprocessor [6] or equivalent.
+
+$$
+\begin{array}{c}{{F_{R d t n}\left(t\right)=-\displaystyle\int_{}^{t}{{K_{R d t n}}\left(t-\tau\right){{\dot{q}}_{P t f m}}\left(\tau\right)d\tau}}}\\ {{\dot{x}_{R d t n}={A_{R d t n}}x_{R d t n}+{B_{R d t n}}{{\dot{q}}_{P t f m}}}}\\ {{F_{R d t n}\cong{C_{R d t n}}x_{R d t n}}}\end{array}
+$$
+
+The linear state-space-based wave-excitation model was newly introduced in HydroDyn within OpenFAST with the addition of the linearization functionality for FOWTs. Similar to the linear state-space-based wave-radiation model, the model approximates the infinite integral shown in Eq. (3) with a linear state-space model as shown in Eq. (4), where $F_{E x c t n}$ is the wave-excitation loads under unidirectional waves; $K_{E x c t n}$ is the incident wave-excitation kernels; $\zeta$ is the wave elevation; $\zeta_{c}$ is a time-shifted wave elevation; $x_{E x c t n}$ is the waveexcitation states; and $A_{E x c t n}$ , $B_{E x c t n}$ , and $C_{E x c t n}$ are the matrices of the linear state-space wave-excitation model. Further details on this new linear state-space-based waveexcitation model are provided in Annex A. Within HydroDyn, the user can now select either a Fourier-transform-based waveexcitation model, which is implemented with discrete Fourier transforms (which is not conducive to linearization), or the new linear state-space-based wave-excitation model.
+
+$$
+\begin{array}{c}{{\displaystyle F_{E x c t n}\left(t\right)=\int_{\phantom{\bigg|}E_{E x c t n}}\left(t-\tau\right)\zeta\left(\tau\right)d\tau}}\\ {{\displaystyle\dot{x}_{E x c t n}=A_{E x c t n}x_{E x c t n}+B_{E x c t n}\zeta_{c}\vphantom{\bigg|}}}\\ {{\displaystyle F_{E x c t n}\cong C_{E x c t n}x_{E x c t n}}}\end{array}
+$$
+
+When linear state-space-based wave excitation or wave radiation are enabled in the potential-flow solution, continuoustime states associated with the excitation and radiation damping are included in the linearized HydroDyn model. Linearization is permitted only in still water $\left(\left.\zeta_{c}\right|_{o p}=x_{E x c t n}\right|_{o p}=F_{E x c t n}\right|_{o p}=0\:)^{3}$ and is not permitted when
+
+Fourier-transform-based wave-excitation, convolution-based wave-radiation, wave directional spreading, second-order wave kinematics, or second-order potential-flow loads are enabled. The linearized form of the HydroDyn state and output equations is given by Eq. (5), where ∆x(HD)= $\varDelta x^{(H D)}=\left\{{\varDelta x_{E x c t n}}\atop{\varDelta x_{R d t n}}\right\},$ $\varDelta u^{(H D)}$ includes contributions from $\varDelta\zeta$ and $\varDelta\dot{q}_{P t f m}$ , and ∆y(HD) includes contributions from ∆FExctn and ∆FRdtn.
+
+$$
+\begin{array}{r c l}{{}}&{{\Delta\dot{x}^{(H D)}=A^{(H D)}A x^{(H D)}+B^{(H D)}A u^{(H D)}}}\\ {{}}&{{}}&{{A y^{(H D)}=C^{(H D)}A x^{(H D)}+D^{(H D)}A u^{(H D)}}}\end{array}
+$$
+
+The continuous-state matrix, input matrix, continuous-state output matrix, and the input-transmission matrix are the Jacobians of the state and output equations relative to the states and inputs about the OP, computed analytically for the Jacobians of the state equations and numerically via a centraldifference perturbation technique for the Jacobians of the output equations as shown in Eq. (6), where $X^{(H D)}$ is the continuous-state functions and $Y^{(H D)}$ is the output functions of HydroDyn. For inputs that are rotations in 3D (i.e., $u=A$ and $\varDelta u=\varDelta\vec{\theta}$ ), it is implied that $u\big|_{o p}^{\left(H D\right)}+\varDelta u^{\left(H D\right)}$ uop , $u\big|_{o p}^{\left(H D\right)}-\varDelta u^{\left(H D\right)}$ , and $2{\varDelta}u^{\left(H D\right)}\;\;\;\;\;\;\mathrm{are}\;\;\;\;\;\mathrm{written}\;\;\;\;\;\mathrm{as}\;\;\;\;\;\;{\varLambda}_{\left\vert{\vphantom{\Biggl(}}^{H D\right)}\right.}^{\left(H D\right)}\;f_{A A}\left(A{\vec{\theta}}^{\left(H D\right)}\right),$ $A\big|_{o p}^{(H D)}\,f_{A A}\big(-\varDelta\vec{\theta}^{(H D)}\big)$ , and $2\varDelta\vec{\theta}^{(H D)}$ , respectively. In Eq. (6b), the column order is dictated by the order of HydroDyn inputs, and $\boldsymbol{\theta}$ are appropriately sized zero matrices. For the numerically computed Jacobians, the default perturbation sizes are hard-coded within HydroDyn (but can be customized by recompiling) and are ${\approx}0.035\mathrm{~m~}$ for translational states, $2^{\circ}$ for rotational states and inputs, and fractions of the water depth for translational inputs. The Jacobians contain the linearized contributions from 1) state-space-based wave excitation, 2) hydrodynamic added mass, 3) state-space-based wave-radiation damping, 4) hydrostatic restoring, and 5) linearized viscous drag. The latter is found by numerically linearizing the viscousdrag term from the relative form of Morison’s equation in the strip-theory solution (plus member-end effects) as in Eq. (6), not based on a stochastic linearization method (e.g., [7]).
+
+The $\mathrm{MAP++}$ (MAP) module has constraint (algebraic) states, there are no restrictions to linearization, and the linearized form is given by $\varDelta y^{\left(M A P\right)}=D^{\left(M A P\right)}\varDelta u^{\left(M A P\right)}$ . The input-transmission matrix is the Jacobian of the output equations relative to the inputs about the OP, computed numerically via a central-difference perturbation technique as shown in Eq. (7), where $Y^{(M A P)}$ is the output functions of $\mathrm{MAP++}$ . For inputs that are rotations in 3D (i.e., $u=A$ and $\begin{array}{r l r}&{\begin{array}{r l r}{\boldsymbol{\varDelta}u=\boldsymbol{\varDelta}\vec{\theta}\,\,\mathrm{)},\quad\mathrm{~it}\quad\mathrm{~is~}\quad\mathrm{~implied~}\quad\mathrm{~that~}\quad u\bigl|_{o p}^{\big(M A P\big)}+\varDelta u^{\big(M A P\big)}\,,}\end{array}}\\ &{\begin{array}{r l r}{u_{o p}^{\big(M A P\big)}-\varDelta u^{\big(M A P\big)}\,,\quad}&{\mathrm{and}\quad}&{\mathrm{2}\varDelta u^{\big(M A P\big)}\quad}&{\mathrm{are~}\quad\mathrm{~written~}\quad\mathrm{~as~}}\end{array}}\\ &{\begin{array}{r l r}{\boldsymbol{\varLambda}_{o p}^{\big(M A P\big)}\,f_{\_{\Delta\boldsymbol{A}}}\left(\boldsymbol{\varDelta}\vec{\theta}^{\big(M A P\big)}\right),\quad}&{\boldsymbol{\varLambda}_{\big|o p}^{\big(M A P\big)}\,f_{\_{\Delta\boldsymbol{A}}}\left(-\varDelta\vec{\theta}^{\big(M A P\big)}\right),\quad}&{\mathrm{and}}\end{array}}\end{array}$ shown in Eqs. (6) and (7) of [2]. This difference in implementation was needed in $\mathrm{MAP++}$ because the residual of constraint equations, $Z^{(M A P)}$ , is not readily available internally within MAP $^{++}$ . The default perturbation sizes are hard-coded within $\mathrm{MAP++}-$ but can be customized by recompiling—and are fractions of the water depth for the inputs. To obtain an accurate Jacobian, it is important that the perturbation sizes are not too large to make a taut mooring line go slack or a slack line go taut, but not too small to prevent the tolerance used to solve the implicit constraint equations from effecting the numerically computed derivative.
+
+$2\varDelta\vec{\theta}^{(M A P)}$ , respectively. The Jacobian contains the linearized contribution from mooring restoring.
+
+# MODULE-TO-MODULE, INPUT-OUTPUT COUPLING RELATIONSHIPS LINEARIZATION
+
+$$
+\begin{array}{r l}&{\mathcal{A}^{(l t D)}=\frac{\hat{\mathcal{A}}\mathcal{X}}{\hat{\mathcal{Z}}\mathcal{X}}\Bigg|_{\varphi=\ D}^{(l t D)}=\Bigg[\frac{\mathcal{A}_{\hat{X}\kappa\sin\varphi}}{\mathcal{A}_{\hat{X}\kappa\sin\varphi}}\Bigg]}\\ &{\mathcal{B}^{(l t D)}=\frac{\hat{\mathcal{B}}\mathcal{X}}{\hat{\mathcal{C}}\mathcal{M}}\Bigg|_{\varphi=\ D}^{(l t D)}=\Bigg[\theta\quad\theta\quad\theta\quad\theta\quad\beta_{E\kappa\omega}\Bigg]}\\ &{\mathcal{C}^{(l t D)}=\frac{\hat{\mathcal{A}}\mathcal{Y}}{\hat{\mathcal{Z}}\mathcal{X}}\Bigg|_{\varphi=\ }^{(l t D)}=\frac{\displaystyle{\cal Y}\left(x\right)_{\varphi}+\varDelta x\,u_{[l]_{\varphi}}\,,I[\mathrm{i}_{\varphi}]}{24\Delta\pi}-\gamma\left(x\right)_{\varphi}-\varDelta x\,u_{[l_{\varphi}}\,,I_{\log}^{\prime}\right)\Bigg|^{(l t D)}}\\ &{\mathcal{D}^{(l t D)}=\frac{\hat{\mathcal{B}}\mathcal{Y}}{\hat{\mathcal{C}}\mathcal{M}}\Bigg|_{\varphi=\varphi}^{(l t D)}=\frac{\displaystyle{\cal Y}\left(x\right)_{\varphi}\,,\,u_{[l_{\varphi}}\,,\beta_{\textup{d}^{\prime}}\,+\,\mathrm{d}u_{t},I_{\varphi}]}{24\Delta\mu}\gamma\left(x\right)_{\varphi}\cdot u_{[l_{\varphi}}^{\prime}\,,-\,\mathrm{d}u_{t},I_{\varphi]}\Bigg|^{(l t D)}}\end{array}
+$$
+
+Most module inputs and outputs in OpenFAST reside on spatial boundaries, which in the FAST modularization framework are defined in terms of a mesh. The module-tomodule, input-output coupling relationships in the OpenFAST glue code are algebraic, and include spatial mesh-to-mesh mapping. The general approach to linearizing the mesh mapping within the module-to-module, input-output coupling relationships within OpenFAST, including rotations, is detailed in [2].
+
+The linearized input-output transformation functions, $U$ , are given by Eq. (8) from [2] or Eq. (15) from [4], repeated in Eq. (8) for convenience.
+
+$$
+D^{(M P)}=\frac{d Y}{d u}\Bigg|_{o p}^{\left(M A P\right)}=\frac{Y\left(z\big|_{o p}+{\varDelta z},u\big|_{o p}+{\varDelta u},t\big|_{o p}\right)-Y\left(z\big|_{o p}-{\varDelta z},u\big|_{o p}-{\varDelta u},t\big|_{o p}\right)}{2{\varDelta u}}\Bigg|^{\left(M P\right)}
+$$
+
+While $\mathrm{MAP++}$ and AeroDyn are both modules with constraint states, their numerical linearization approach differs. Equation (7) involves the total derivative of the $\mathrm{MAP++}$ output equations relative to the inputs about the OP, where it is implied that the constraint perturbations, ∆z(MAP) and $-\varDelta z^{(M A P)}$ , are tied to the input perturbations, $\varDelta u^{(M A P)}$ and $-\varDelta u^{\left(M A P\right)}$ , by solving the implicit constraint equations, $0=Z\Bigl(z\big|_{o p}+{\cal A}z,u\big|_{o p}+{\cal A}u,t\big|_{o p}\Bigr)\Bigr|^{(M A P)}$ and $0=Z\Bigl(z\big|_{o p}-varDelta z,u\big|_{o p}-\varDelta u,t\big|_{o p}\Bigr)\Bigr|^{(M A P)}$ , respectively, essentially internally eliminating the constraint states within $\mathrm{MAP++}$ during the linearization process. This is in contrast to the approach taken to linearize AeroDyn, whereby the Jacobians $\left.\frac{\partial Z}{\partial z}\right|_{o p}^{\left(A D\right)},\;\;\left.\frac{\partial Z}{\partial u}\right|_{o p}^{\left(A D\right)},\;\;\left.\frac{\partial Y}{\partial z}\right|_{o p}^{\left(A D\right)},\;\;\mathrm{and}\;\;\left.\frac{\partial Y}{\partial u}\right|_{o p}^{\left(A D\right)}$ are computed separately and algebraically manipulated to compute D as
+
+$$
+{\cal{O}}=\frac{\partial U}{\partial\tilde{u}}\Bigg\vert_{o p}\,\varDelta u+\frac{\partial U}{\partial y}\Bigg\vert_{o p}\,\varDelta y\,\,\mathrm{with}\,\left\vert\frac{\partial U}{\partial\tilde{u}}\right\vert_{o p}\Bigg\vert\neq0
+$$
+
+As is evident from Table 1, the InflowWind, ServoDyn, ElastoDyn, AeroDyn, BeamDyn, HydroDyn, and $\mathrm{MAP++}$ modules were developed so that for the most part—other than mapping between independent spatial discretizations—the input of one module equals the output of another. It follows
+
+that
+
+$$
+U=\left\{U^{\left(J W\right)}\right\},\qquad\quad A u=\left[\begin{array}{c}{{\left(U^{\left(J W\right)}\right)}}\\ {{U^{\left(S\times J D\right)}}}\\ {{U^{\left(E D\right)}}}\\ {{U^{\left(B D\right)}}}\\ {{\left(U^{\left(J D\right)}\right)}}\\ {{U^{\left(4D\right)}}}\\ {{U^{\left(H D\right)}}}\\ {{U^{\left(M D\right)}}}\\ {{U^{\left(M A P\right)}}}\end{array}\right],\qquad\begin{array}{c}{{\left[A u^{\left(J W\right)}\right]}}\\ {{\left(A u^{\left(S r D\right)}\right)}}\\ {{A u^{\left(B D\right)}}}\\ {{A u^{\left(B D\right)}}}\\ {{\left(A u^{\left(A D\right)}\right)}}\\ {{\left(A u^{\left(A D\right)}\right)}}\\ {{\left(A u^{\left(A D\right)}\right)}}\end{array}\right],
+$$
+
+
+
+from Eq. (8) for these seven modules are given by Eq. (9), where $I$ and $\boldsymbol{\theta}$ are appropriately sized identity and zero matrices and the sub-Jacobian matrices are composed of $I\mathrm{~s~}$ , $\theta\;\mathbf{s}.$ , and the linearized matrices from the mapping transfers given in [2].
+
+
+
+Because of their large size, the sub-Jacobian matrices are not shown here, but are described qualitatively instead:
+
+The first, second, fourth, and fifth equations of Eq. (8) for InflowWind inputs, ${\boldsymbol{\theta}}={\boldsymbol{U}}^{(I J W)}$ , ServoDyn inputs, $\boldsymbol{\theta}=\boldsymbol{U}^{(S r v D)}$ , BeamDyn inputs, $\boldsymbol{{\cal O}}=\boldsymbol{{U}}^{(B D)}$ , and AeroDyn inputs, $0=U^{(A D)}$ , are described in [2,3]. The third equation of Eq. (8) for ElastoDyn inputs, $0=U^{(E D)}$ , are also described in [2,3], except that the new floating offshore linearization functionality now also includes terms expressing the contributions of HydroDyn and $\mathrm{MAP++}$ . That is, this equation expresses that the applied point force and moment perturbations distributed along the blades and tower as input to ElastoDyn are derived from the aerodynamicapplied line (per-unit length) force and moment perturbations distributed along the blades and tower as output from AeroDyn. This linearized load-mapping transfer also depends on the translational-displacement perturbations of analysis nodes along the blades and tower output from ElastoDyn. And when BeamDyn is enabled, point force and moment perturbations on the hub as input to ElastoDyn are derived from the bladeroot reaction point force and moment output from BeamDyn, which also depends on the translationaldisplacement perturbation of the hub reference point output from ElastoDyn. Additionally, the blade-pitchangle-commands, nacelle-yaw-moment, and generator-torque perturbations as input to ElastoDyn are derived from the equivalent outputs from ServoDyn. For the floating offshore functionality, point force and moment perturbations on the platform as input to ElastoDyn are derived from the hydrodynamic-applied line (per-unit length) and point force and moment perturbations distributed along the floating platform as output from HydroDyn and the reaction point forces (tensions) lumped at each fairlead as output from $\mathrm{MAP++}$ , which also depends on the translational-displacement perturbation of the platform reference point output from ElastoDyn.
+
+The sixth equation of Eq. (8) for HydroDyn inputs, $O=U^{(H D)}$ expresses that the translational displacement, orientation, translational and rotational velocity, and translational and rotational acceleration perturbations of analysis nodes distributed along the floating platform as input to HydroDyn are derived from motion outputs from ElastoDyn at the platform reference point. The seventh equation of Eq. (8) for $\mathrm{MAP++}$ inputs, $\boldsymbol{\theta}=\boldsymbol{U}^{(M A P)}$ expresses that the translational displacement perturbations of each fairlead as input to $\mathrm{MAP++}$ are derived from the motion outputs from ElastoDyn at the platform reference point.
+
+The Jacobian, $\left.\frac{\partial{\cal U}}{\partial\tilde{u}}\right|_{o p}$ , has ones along its entire diagonal and it is easily shown that its determinant from Eq. (8) is nonzero.
+
+# FINAL MATRIX ASSEMBLY
+
+Once all individual modules and input-output relationships are linearized about the OP, the linearized model of the complete coupled system can be assembled. Linearization of the full-system model produces a linear state-space model representation of the complete nonlinear system about the OP, including the influence of system state and input perturbations on the system response and outputs. The general linearized form of the complete coupled system is given by Eqs. (18) and form for the OpenFAST features linearized to date (without discrete-time states and with ElastoDyn, BeamDyn, and HydroDyn as the only modules with continuous-time states) is given by Eq. (10), where $\boldsymbol{\varDelta u}^{+}$ is the additional input perturbations (explained further in [4]).
+
+$$
+\begin{array}{r}{\varDelta\dot{x}=A A x+B\varDelta u^{+}}\\ {\varDelta y=C\varDelta x+D\varDelta u^{+}}\end{array}
+$$
+
+The full-system state-space matrices are given in Eq. (11), where $G\bigr|_{o p}$ —explained more in [4] for the general case—for the OpenFAST linearization to date is given by Eq. (12). The matrix $G\Big\vert_{o p}$ has ones along its entire diagonal and it is easily shown that its determinant from Eq. (12) is nonzero, which means that the matrix inverse, $\left[G\big|_{o p}\right]^{-I}$ from Eq. (11), exists and is bounded in the neighborhood around the OP.
+
+The input-transmission matrices impact all matrices of the linearized coupled system, highlighting the important role played by direct feedthrough of input to output in the coupled system response. For example, while the continuous-state matrix of ElastoDyn, $\boldsymbol{A}^{(E D)}$ , contains mass, stiffness, and damping only directly associated with the structural model, the full-system continuous state matrix, $A$ , contains mass, stiffness, and damping associated with coupled aero-hydroservo-elastics, including hydrodynamic added mass, waveradiation and viscous damping, and hydrostatic and mooring restoring for FOWTs.
+
+When the linearized full-system matrices $A\,,\,\,B\,,\,\,C$ , and $D$ are exported to a file by OpenFAST, the additional input perturbations, ${\varDelta}{u}^{+}$ , can be chosen by the user to be: 1) the inputs of all modules, 2) none of the module inputs (removing $B$ and $D$ from the file), or 3) a standard subset of these inputs, which include the standard wind turbine control inputs of nacelle-yaw moment, generator torque, and blade-pitchangle commands (both independent and rotor-collective); the standard wind-inflow disturbances of horizontal wind speed, power-law shear exponent, and wind-propagation direction; and the standard incident-wave disturbance of wave elevation. Likewise, the output perturbations, $_{\varDelta y}$ , can be chosen by the user to be: 1) the outputs of all modules, 2) none of the module outputs (removing $C$ and $D$ from the file), or 3) only the subset of output variables selected by the user through the OpenFAST module input file(s). Regardless of what the user selects to be exported to a file, all of the module inputs and outputs are used to form the linearized full-system matrices in Eq. (11), but only a subset of these matrices are exported based on the user selection.
+
+$$
+{\left[\begin{array}{l l l l l l l l}{\left|{\begin{array}{l l l l l l l l l}{\alpha}&&{\beta^{\operatorname{aff}}}&&&{0}&&{0}&{\ddots}&{0}\\ {\alpha}&{\beta^{\operatorname{aff}}}&{0}&&{0}&&{\beta^{\operatorname{aff}}\left|{\vphantom{\begin{array}{l l l l l l l l}{0}}&&{0}&&{0}&&{0}\\ {0}&&{\gamma^{\operatorname{aff}}}&{0}&&{0}&&{0}&{0}\\ {0}&{0}&&{\gamma^{\operatorname{aff}}}&{\left|{\begin{array}{l l l l l l l l}{\alpha}&{0}&&{\beta^{\operatorname{aff}}}&&{0}&&{0}\\ {0}&&{0}&&{0}&&{\theta^{\operatorname{aff}}}&{\left|{\vphantom{\begin{array}{l l l l l l l l}{{\ddots}}&{0}&&{0}&{0}\\ {0}&&{0}&&{0}&{0}&{\ddots}&&{0}\\ {0}&&{0}&&{0}&&{\gamma^{\operatorname{aff}}}&{\left|{\vphantom{\begin{array}{l l l l l l l l l}{0}}&&{0}&&{0}&{0}\\ {0}&&{0}&&{\gamma^{\operatorname{aff}}}&&{\left|{\begin{array}{l l l l l l l l l}{0}&&{0}&&{0}&{0}\\ {0}&&{0}&&{0}&{\ddots}&&{0}\\ {0}&&{0}&&{0}&{0}&&{0}&{\gamma^{\operatorname{aff}}}&{\left|{\vphantom{\begin{array}{l l l l l l l l l}{0}&{0}&{0}&{0}&{0}\\ {0}&&{0}&{0}&{\ddots}&&{0}\\ {0}&&{0}&{0}&&{0}&{0}&{\gamma^{\operatorname{aff}}}&{\left|{\vphantom{\begin{array}{l l l l l l l l l}{0}&{0}&{0}&{0}&{0}\\ {0}&&{0}&{0}&{0}&{\gamma^{\operatorname{aff}}}&&{\left|{\begin{array}{l l l l l l l l l}{0}&&{0}&{0}&{0}\\ {0}&&{0}&{0}&{0}&{0}\\ {0}&&{0}&{\gamma^{\operatorname{aff}}}&{0}&&{0}&{0}\\ {0}&{0}&&{0}&{0}&{
+$$
+
+
+
+# CONCLUSIONS
+
+Linearization of the underlying nonlinear wind-system equations is often important for understanding the system response and exploiting well-established methods and tools for analyzing linear systems. This paper has presented the development of the new linearization functionality of the opensource engineering tool OpenFAST for FOWTs, as well as the concepts and mathematical background needed to understand and apply it. Details have been presented on finding an OP, linearizing the underlying nonlinear equations of each module about the OP, linearizing the module-to-module input-output coupling relationships about the OP, and combining all linearized matrices into the full-system linear state-space model. The new linearization functionality enables the contributions from state-space-based wave excitation, hydrodynamic added mass, state-space-based wave-radiation damping, hydrostatic restoring, linearized viscous drag, and linearized mooring restoring for FOWTs to be included in the full linearized system.
+
+Unfortunately, the implementation at the time of this writing has not yet been completed enough to produce results. Results will be presented in future work to highlight the functionality and verify the implementation (e.g., for a sample case whereby the natural frequencies and damping of the NREL 5-MW baseline wind turbine atop the OC3-Hywind spar buoy will be calculated as a function of rotational speed and presented in Campbell-diagram form).
+
+Beyond the developments implemented here, potential future OpenFAST linearization enhancements include: 1) introducing routines to find a static equilibrium, steady-state, or periodic steady-state condition for improved OP determination; 2) enabling the linearization of the tuned-mass-damper (TMD)
+
+models and external user-specified controllers within ServoDyn for improved controls development and analysis; 3) implementing a form of unsteady airfoil aerodynamics, including dynamic stall, within AeroDyn amendable to linearization (e.g., from [8]) and including that model within the linearization of AeroDyn for advanced aerodynamic and stability analysis; 4) linearizing the offshore functionality of OpenFAST for bottom-fixed offshore wind systems; and 5) linearizing new features as the OpenFAST models are further developed in the future going forward.
+
+# ACKNOWLEDGMENTS
+
+The authors would like to thank Bonnie Jonkman of Envision Energy USA, Ltd. for her contributions to the development of the FAST modularization framework and to implementation of the land-based linearization functionality.
+
+This work was authored by the National Renewable Energy Laboratory, operated by the Alliance for Sustainable Energy, LLC, for the U. S. Department of Energy (DOE) under Contract No. DE-AC36-08GO28308. Funding was provided by the U.S. Department of Energy Office of Energy Efficiency and Renewable Energy. The views expressed in the article do not necessarily represent the views of the U.S. Department of Energy or the U.S. government. The U.S. government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this work, or allow others to do so, for U.S. government purposes.
+
+# REFERENCES
+
+[1] Web page https://nwtc.nrel.gov/OpenFAST (accessed May 7, 2018).
+
+[2] Jonkman, J. M. and Jonkman, B. J. “FAST modularization framework for wind turbine simulation: full-system linearization.” Journal of Physics: Conference Series, The Science of Making Torque from Wind (TORQUE 2016), 5–7 October 2016, Munich, Germany [online journal]. 082010. Vol. 753, 2016. URL: http://iopscience.iop.org/article/10.1088/1742- 6596/753/8/082010/pdf; NREL/CP-5000-67015. Golden, CO: National Renewable Energy Laboratory.
+[3] Jonkman, J. M.; Jonkman, B. J.; and Platt, A. “FullSystem Linearization for Wind Turbines with Aeroelastically Tailored Rotor Blades in OpenFAST.” (in preparation).
+[4] Jonkman, J. M. “The New Modularization Framework for the FAST Wind Turbine CAE Tool.” $5l^{s t}$ AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, 7–10 January 2013, Grapevine (Dallas/Ft. Worth Region), TX [online proceedings]. URL: http://arc.aiaa.org/doi/pdf/10.2514/6.2013-202. AIAA-2013-0202. Reston, VA: American Institute of Aeronautics and Astronautics, January 2013; NREL/CP5000-57228. Golden, CO: National Renewable Energy Laboratory.
+[5] Duarte, T.; Alves, M.; Jonkman, J.; and Sarmento, A. “State-Space Realization of the Wave-Radiation Force within FAST.” $32^{n d}$ International Conference on Ocean, Offshore, and Arctic Engineering (OMAE2013), 9–14 June 2013, Nantes, France [DVD-ROM]. OMAE2013- 10375. Houston, TX: The American Society of Mechanical Engineers (ASME International) Ocean, Offshore and Arctic Engineering (OOAE) Division, June 2013; ISBN 978-0-7918-5542-3; NREL/CP-5000-58099. Golden, CO: National Renewable Energy Laboratory.
+[6] Web page https://nwtc.nrel.gov/SS_Fitting (accessed May 7, 2018).
+[7] Borgman, L. E. “Spectral Analysis of Ocean Wave Forces on Piling.” Journal of the Waterways and Harbors Division, Vol. 93, Issue 2, pp. 129-156, 1967.
+[8] Hansen, M. H.; Guanaa, M.; and Madsen, H. A. A Beddoes-Leishman type dynamic stall model in statespace and indicial formulations. Risø-R-1354 (EN). Roskilde, Denmark: Risø National Laboratory. 2004.
+
+# ANNEX A
+
+# LINEAR STATE-SPACE-BASED WAVE-EXCITATION MODEL
+
+In the frequency domain, the potential-flow-based waveexcitation loads (three forces and three moments) can be expressed for unidirectional waves by Eq. (A.1), where $X\big(\omega,\beta\big)$ is the frequency- $\left(\,\omega\,\right)$ and direction- $\left(\,\beta\,\right)$ dependent wave-excitation loads on the floating platform normalized per unit wave amplitude (derived, e.g., from a potential-flow-based wave-body interaction solver WAMIT) and $\zeta(\omega)$ is the waveelevation spectrum for waves propagating in direction $\beta$ (both are complex-valued to include phases). For time-domain simulations, the spectral wave data and wave direction are model inputs and the phases of the various wave components are selected randomly. The time-domain wave-excitation loads, $F_{E x c t n}\left(t\right)$ , are usually precomputed at model initialization by calculating the inverse Fourier transform shown in Eq. (A.2), where $j=\sqrt{-I}$ is the imaginary number (HydroDyn uses a discrete Fourier transform in place of the continuous-time transform shown).
+
+$$
+F_{E x c t n}\left(\omega\right)=X\big(\omega,\beta\big)\zeta\big(\omega\big)
+$$
+
+$$
+F_{E x c t n}\left(t\right)=\frac{I}{2\pi}\int_{-\infty}^{\infty}X\big(\omega,\beta\big)\zeta\left(\omega\right)\!e^{j\omega t}d\omega
+$$
+
+Equivalently, as shown in [A.1], the wave-excitation load components can be calculated through the infinite convolution integral given by Eq. (3), where the incident wave-excitation kernels derived from Eq. (A.3), represent the incident waveexcitation impulse-response functions (IRF) that depend only on platform parameters. The goal is to use a linear state-space equation to approximate Eq. (3), with the time-domain wave elevation, $\zeta\!\left(t\right)$ , acting as the input, and the time-domain wave-excitation loads, $F_{E x c t n}\left(t\right)$ , as the output.
+
+$$
+K_{E x c t n}\left(t\right)=\frac{I}{2\pi}\int_{-\infty}^{\infty}X\big(\omega,\beta\big)e^{j\omega t}d\omega
+$$
+
+Fitting a linear state-space model to the infinite convolution integral is complicated by the noncausality of this system. This issue has been studied extensively in [A.2] and [A.3], where a time-shifting method is used to arrive at a causal transfer function between wave elevation and wave-excitation loads. We follow the same method used in [A.2] and [A.3] here. Figure A.1 shows an example platform-pitch response (blue curve) resulting from a wave-elevation impulse at time $t=\partial$ s. Clearly, the system is noncausal because its response is nonzero for negative time even though the impulse acts at $t=\partial$ s. The system response for negative time is as large as the response for positive time, so it cannot be ignored.
+
+
+FIGURE A.1. NONCAUSAL AND CAUSAL IMPULSERESPONSE FUNCTIONS.
+
+An approximately causal IRF is obtained through a timeshifting method. Figure A.1 shows that for all times less than $-t_{c}$ ( $t_{c}\approx I2$ s for this example), the IRF is approximately zero. A time delay, $t_{c}$ , is introduced to shift the original IRF to obtain a “causalised” IRF, $K_{E x c t n_{c}}\left(t\right)=K_{E x c t n_{c}}\left(t-t_{c}\right).$ . Figure A.1 shows that this time-shifted IRF (red) is approximately zero for negative time. Equation (A.4), in terms of $K_{E x c t n_{c}}\left(t\right)$ and $\zeta_{c}\left(t\right)$ , follows directly from Eq. (3) and the time-shift properties of Fourier transforms, where $\zeta_{c}\left(t\right)=\zeta\left(t+t_{c}\right)$ is the wave elevation predicted at time $t_{c}$ into the future. To implement this approach, a wave-prediction model must be used to predict the wave elevation for future times. Within OpenFAST simulations, this is not an issue because the wave elevation is precalculated at model initialization, so the future wave elevations are already available. In practice, the future wave evaluation could be estimated (e.g., by measuring the wave elevation with a buoy and approximating the wave evolution).
+
+$$
+F_{E x c t n}\left(t\right)=\int_{-\infty}^{\infty}K_{E x c t n_{c}}\left(t-\tau\right)\zeta_{c}\left(\tau\right)d\tau
+$$
+
+Now the goal is to fit a linear state-space model to the causalised IRF, $K_{E x c t n_{c}}\left(t\right)$ , based on input $\zeta_{c}\left(t\right).$ . For a statespace system of the form shown in Eq. (4), it can be shown that the IRF is $C_{E x c t n}\ e^{A_{E x c t n}t}\ B_{E x c t n}$ , where $e^{\b{A}_{E x c t n}t}$ is the matrixexponential function [A.4]. The goal is to determine $A_{E x c t n}$ , $B_{E x c t n}$ , and $C_{E x c t n}$ so that $C_{E x c t n}\ e^{A_{E x c t n}t}\ B_{E x c t n}$ approximates $K_{E x c t n_{c}}\left(t\right)$ in terms of some error norm. This is usually done by minimizing a cost function, such as given by Eq. (A.5), where $G\big(t_{k}\big)$ is a weighting function to be chosen, evaluated at specific sampling times, $t_{k}$ , for $k=I,2,\dotsc,m$ . At least $m$ discrete values of KExctn a re assumed to be known either from theory or experiment.
+
+$$
+Q=\sum_{k=I}^{m}G\!\left(t_{k}\right)\!\left[K_{E x c t n_{c}}\left(t_{k}\right)\!-C_{E x c t n}~e^{A_{E x c t n}t}~B_{E x c t n}\right]^{2}
+$$
+
+Several methods exist for performing this systemidentification step. For this work, we rely on the use of a linear state-space realization of the impulse response using Hankel singular value decomposition and MATLAB’s SystemIdentification-Toolbox command imp2ss [A.5]. The state-space model structure is selected so that the infinite-frequency limit of the identified model in the frequency domain is zero, forcing the output transmission $(\,D\,)$ matrix in the general state-space description to be equal to the zero matrix [A.5]. A balanced model-truncation method is used to arrive at a low-order statespace approximation for the wave-excitation loads.
+
+This general method for generating a linear state-space model of the wave excitation for inclusion into OpenFAST is being developed and programmed into a MATLAB script. Further details, with example cases and a user’s guide for this script will be described in a future publication.
+
+[A.1] Jonkman, J. M. Dynamics Modeling and Loads Analysis of an Offshore Floating Wind Turbine. Ph.D. Thesis. Department of Aerospace Engineering Sciences, University of Colorado, Boulder, CO, 2007; NREL/TP500-41958. Golden, CO: National Renewable Energy Laboratory.
+[A.2] Yu, Z. and Falnes, J. “State-space Modeling of a Vertical Cylinder in Heave.” Applied Ocean Research, Vol. 17, Issue 5, pp. 265-275, 1995.
+[A.3] Lemmer, F., Raach, S., Schlipf, D., and Cheng, P. “Parametric Wave Excitation Model for Floating Wind Turbines.” Energy Procedia, Vol. 94, 2016, pp. 290-305; $I3^{t h}$ Deep Sea Offshore Wind R&D Conference, EERA DeepWind’2016, 20-22 January 2016, Trondheim, Norway; DOI: 10.1016/j.egypro.2016.09.186.
+[A.4] Ogata, K. Modern Control Engineering, Englewood Cliffs, NJ: Prentice-Hall Inc., 1990.
+[A.5] Ljung, L. and Singh, R. “Version 8 of the MATLAB System Identification Toolbox.” $I6^{t h}$ IFAC Symposium on System Identification, Brussels, Belgium. ISBN 9783902823069; pp. 1826-1832, 2012; doi:10.3182/20120711-3-BE-2027.00061.
\ No newline at end of file
diff --git a/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/Jonkman-Full-system linearization for floating_origin.pdf b/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/Jonkman-Full-system linearization for floating_origin.pdf
new file mode 100644
index 0000000..75a799a
Binary files /dev/null and b/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/Jonkman-Full-system linearization for floating_origin.pdf differ
diff --git a/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/images/228dd88a114ae163dad8ffb4be4a2a74d82c5bee4a868bbce9f447707580824f.jpg b/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/images/228dd88a114ae163dad8ffb4be4a2a74d82c5bee4a868bbce9f447707580824f.jpg
new file mode 100644
index 0000000..5264513
Binary files /dev/null and b/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/images/228dd88a114ae163dad8ffb4be4a2a74d82c5bee4a868bbce9f447707580824f.jpg differ
diff --git a/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/images/22dc794cd6fdc6791a3412b34d67e99749275b20bb63eee47ea19d5908a94158.jpg b/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/images/22dc794cd6fdc6791a3412b34d67e99749275b20bb63eee47ea19d5908a94158.jpg
new file mode 100644
index 0000000..e5a4948
Binary files /dev/null and b/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/images/22dc794cd6fdc6791a3412b34d67e99749275b20bb63eee47ea19d5908a94158.jpg differ
diff --git a/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/images/63a387a5c4d629058ef1038d42302921b057bc3c4140965b126080002098ca7c.jpg b/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/images/63a387a5c4d629058ef1038d42302921b057bc3c4140965b126080002098ca7c.jpg
new file mode 100644
index 0000000..e1733be
Binary files /dev/null and b/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/images/63a387a5c4d629058ef1038d42302921b057bc3c4140965b126080002098ca7c.jpg differ
diff --git a/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/images/a831b7ac45c7ef3e2ffe228fb1164fc895012b24d942fac4f79d60ba618f2485.jpg b/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/images/a831b7ac45c7ef3e2ffe228fb1164fc895012b24d942fac4f79d60ba618f2485.jpg
new file mode 100644
index 0000000..393f154
Binary files /dev/null and b/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/images/a831b7ac45c7ef3e2ffe228fb1164fc895012b24d942fac4f79d60ba618f2485.jpg differ
diff --git a/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/images/e3cf5199761175d226cc33794ee9ebb7a89b2f9cf0f39528b2eaeef3d2ad016d.jpg b/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/images/e3cf5199761175d226cc33794ee9ebb7a89b2f9cf0f39528b2eaeef3d2ad016d.jpg
new file mode 100644
index 0000000..00b03a5
Binary files /dev/null and b/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/images/e3cf5199761175d226cc33794ee9ebb7a89b2f9cf0f39528b2eaeef3d2ad016d.jpg differ
diff --git a/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/images/f47cb8194eff840ddeca7a7167915ca9ff6fca5b6c65e89e27039e48f46917cc.jpg b/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/images/f47cb8194eff840ddeca7a7167915ca9ff6fca5b6c65e89e27039e48f46917cc.jpg
new file mode 100644
index 0000000..9c23f64
Binary files /dev/null and b/线性化求解器/参考文献/Jonkman-Full-system linearization for floating/auto/images/f47cb8194eff840ddeca7a7167915ca9ff6fca5b6c65e89e27039e48f46917cc.jpg differ
diff --git a/线性化求解器/参考文献/Skjoldan-Implicit Floquet analysis of wind tur/auto/Skjoldan-Implicit Floquet analysis of wind tur.md b/线性化求解器/参考文献/Skjoldan-Implicit Floquet analysis of wind tur/auto/Skjoldan-Implicit Floquet analysis of wind tur.md
new file mode 100644
index 0000000..9c5c302
--- /dev/null
+++ b/线性化求解器/参考文献/Skjoldan-Implicit Floquet analysis of wind tur/auto/Skjoldan-Implicit Floquet analysis of wind tur.md
@@ -0,0 +1,377 @@
+# Implicit Floquet analysis of wind turbines using tangent matrices of a non-linear aeroelastic code
+
+P. F. Skjoldan1 and M. H. Hansen2
+
+1 Loads, Aerodynamic and Control, Siemens Wind Power A/S, DK-2630 Taastrup, Denmark
+2 Wind Energy Division, National Laboratory for Sustainable Energy, Risø DTU, DK-4000 Roskilde, Denmark
+
+# ABSTRACT
+
+The aeroelastic code BHawC for calculation of the dynamic response of a wind turbine uses a non-linear finite element formulation. Most wind turbine stability tools for calculation of the aeroelastic modes are, however, based on separate linearized models. This paper presents an approach to modal analysis where the linear structural model is extracted directly from BHawC using the tangent system matrices when the turbine is in a steady state. A purely structural modal analysis of the periodic system for an isotropic rotor operating at a stationary steady state was performed by eigenvalue analysis after describing the rotor degrees of freedom in the inertial frame with the Coleman transformation. For general anisotropic systems, implicit Floquet analysis, which is less computationally intensive than classical Floquet analysis, was used to extract the least damped modes. Both methods were applied to a model of a three-bladed $2.3\;\mathrm{MW}$ Siemens wind turbine model. Frequencies matched individually and with a modal identification on time simulations with the non-linear model. The implicit Floquet analysis performed for an anisotropic system in a periodic steady state showed that the response of a single mode contains multiple harmonic components differing in frequency by the rotor speed. Copyright $\copyright$ 2011 John Wiley & Sons, Ltd.
+
+# KEYWORDS
+
+modal analysis; Floquet analysis; rotor dynamics
+
+# Correspondence
+
+P. F. Skjoldan, Loads, Aerodynamic and Control, Siemens Wind Power A/S, Dybendalsvænget 3, DK-2630 Taastrup, Denmark. E-mail: peter.skjoldan@siemens.com
+
+Received 26 June 2010; Revised 7 October 2010; Accepted 12 February 2011
+
+# 1. INTRODUCTION
+
+Today, advanced non-linear finite element codes1–3 are routinely used for load calculations on wind turbines. Most wind turbine stability tools for calculation of the aeroelastic modes are, however, based on separate linearized models. Stability analysis can be divided into three steps: first, a calculation of the steady state; then, a linearization of the equations of motion about the steady state and last, a modal analysis to extract modal frequencies, damping and mode shapes. This paper presents an approach to structural modal analysis applicable to any periodic steady state where the linearization is obtained directly from the non-linear wind turbine aeroelastic code BHawC.3
+
+The equations of motion for a wind turbine operating at a constant mean rotor speed contain periodic coefficients, preventing direct eigenvalue analysis of the system. Most recent wind turbine stability tools $4\mathrm{-}7$ incorporate the Coleman transformation, also known as the multiblade coordinate transformation, which describes the rotor degrees of freedom in the inertial frame. This transformation eliminates the periodic coefficients if the system is isotropic, i.e. the rotor consists of identical symmetrically mounted blades, and the environment conditions are symmetric. Floquet analysis is, however, applicable to anisotropic systems and any periodic steady state. It requires integration of the equations of motion over a period of rotor rotation, as many times as there are state variables in the system. Because of the computational burden of this approach, it has only been applied to reduce or simplify wind turbine models with a limited number of degrees of freedom.8–10 One way to reduce the computation time is to use the Fast Floquet Theory11 where only one third of the integrations are necessary for a three-bladed isotropic rotor. Another way is to use implicit Floquet analysis12 where the least damped modes can be extracted after a limited number of integrations.
+
+Stol et al.13 compare the Floquet analysis with the Coleman transformation approach applied to a periodic steady state, where the remaining periodic coefficients are eliminated by averaging and find small differences in modal frequencies and damping, concluding that it is not necessary to use Floquet analysis.
+
+Another approach to modal analysis is system identification,14–16 which operates on the response from numerical simulations or measurements, and no knowledge of the system equations is needed to extract the modal properties. The accuracy of the methods is, however, limited and depends on the chosen excitation.
+
+In this paper, tangent matrices for mass, damping and stiffness are extracted from the aeroelastic code BHawC. If the system is isotropic and the steady state is stationary, the Coleman transformation is applied before extracting the modal parameters by eigenvalue analysis. For an anisotropic system, implicit Floquet analysis is used for the modal analysis. When the system is isotropic, the response of a single mode contains a single harmonic component for tower degrees of freedom and up to three components for the blades. The response of a single mode in the anisotropic system on both blades and tower contains multiple harmonic components differing in frequency by the rotor speed.
+
+Section 2 of this paper describes the BHawC model, and Section 3 explains the methods for modal analysis, the Coleman transformation approach, the implicit Floquet analysis and also the partial Floquet analysis, a system identification technique. In Section 4, the methods are applied to a model of a wind turbine. Section 5 discusses the approaches, and Section 6 concludes the paper.
+
+# 2. STRUCTURAL MODEL
+
+The BHawC wind turbine aeroelastic code3 is based on a structural finite element model sketched in Figure 1, where the main structural parts, tower, nacelle, shaft, hub and blades, are modelled as two-node 12-degrees of freedom Timoshenko beam elements. The code uses a corotational formulation, where each element has its own coordinate system that rotates with the element. The elastic deformation is described in the element frame, whereas the movement of the element coordinate system accounts for rigid body motion. In this way, a geometrically non-linear model is obtained using linear finite elements.
+
+The configuration of the system, defined by nodal positions $\pmb{p}$ and orientations $\pmb q$ , nodal velocities $\dot{\pmb u}$ (of both positions and orientations) and nodal accelerations $\ddot{u}$ , must satisfy the equilibrium equation given in global coordinates as
+
+$$
+f_{\mathrm{iner}}(\boldsymbol{p},\boldsymbol{q},\dot{\boldsymbol{u}},\ddot{\boldsymbol{u}})+f_{\mathrm{damp}}(\boldsymbol{q},\dot{\boldsymbol{u}})+f_{\mathrm{int}}(\boldsymbol{p},\boldsymbol{q})=f_{\mathrm{ext}}
+$$
+
+where $f_{\mathrm{iner}},f_{\mathrm{damp}},f_{\mathrm{int}}$ and $\pmb{f}_{\mathrm{ext}}$ are the inertial, damping, internal and external force vectors, respectively, and $\dot{\mathrm{O}}=\mathrm{d}/\mathrm{d}t$ denotes a time derivative. The inertial forces depend on the acceleration of the masses, the damping forces are given by viscous damping, the internal forces are due to elastic forces and the external forces contain the aerodynamic forces.17 To find
+
+
+Figure 1. Sketch of the BHawC model substructures.
+
+this equilibrium configuration, increments of the positions and the orientations $\delta\pmb{u}$ , the velocities $\delta\dot{\pmb{u}}$ and the accelerations $\delta\ddot{\pmb{u}}$ are obtained using Newton–Raphson iteration with the tangent relation obtained from the variation of Equation (1) as
+
+$$
+\mathbf{M}(q)\delta{\ddot{u}}+\mathbf{C}(q,{\dot{u}})\delta{\dot{u}}+\mathbf{K}(p,q,{\dot{u}},{\ddot{u}})\delta u=r
+$$
+
+where M, C and $\mathbf{K}$ are the tangent mass, damping/gyroscopic and stiffness matrices, respectively, and $r=f_{\mathrm{ext}}+f_{\mathrm{iner}}+$ $f_{\mathrm{damp}}\!-\!f_{\mathrm{int}}$ is the residual. The stiffness matrix is composed of constitutive, geometric and inertial stiffness. The orientation $\pmb q$ of the nodes is described by quaternions, also known as the Euler parameters,18 a general four-parameter representation equivalent to a triad, which for node number $i$ is updated as
+
+$$
+\pmb q_{i}:=q u a t(\delta\pmb u_{i,\mathrm{rot}})*\pmb q_{i}
+$$
+
+where $\delta\mathbf{\boldsymbol{u}}_{i,\mathrm{rot}}$ contains three rotations that are assumed infinitesimal and thus commute and where this rotation pseudovector is transformed by the function termed quat into a quaternion, which is used to update the nodal quaternion $\pmb q_{i}$ employing the special quaternion product denoted by $^*$ , which maintains the unity of the quaternion. The nodal positions $\pmb{p}$ , the nodal velocities $\dot{\pmb u}$ and the accelerations $\ddot{u}$ are updated by regular addition of the positional part of $\delta\pmb{u},\,\delta\dot{\pmb{u}}$ and $\delta\ddot{\pmb{u}}$ , respectively. All components in $\pmb{p}$ , $\pmb q$ and $\delta\pmb{u}$ are absolute and described in a global frame.
+
+The present work considers small perturbations in position and orientation ${\bf\delta y}$ , velocity $\dot{\mathbf{y}}$ and acceleration $\ddot{\mathbf{y}}$ to a steady state with constant mean rotor speed $\varOmega$ defined by $(p_{\mathrm{ss}},q_{\mathrm{ss}},\dot{\pmb u}_{\mathrm{ss}},\ddot{\pmb u}_{\mathrm{ss}})$ , the steady state positions, orientations, velocities and accelerations, respectively, all periodic with the rotor period $T=2\pi/\varOmega$ . The linearized equations of motion are obtained from equation (2) at $r\approx\theta$ as
+
+$$
+{\bf M}({q}_{\mathrm{ss}})\ddot{\boldsymbol{y}}+{\bf C}({q}_{\mathrm{ss}},\dot{\boldsymbol{u}}_{\mathrm{ss}})\dot{\boldsymbol{y}}+{\bf K}({p}_{\mathrm{ss}},{q}_{\mathrm{ss}},\dot{\boldsymbol{u}}_{\mathrm{ss}},\ddot{\boldsymbol{u}}_{\mathrm{ss}}){\boldsymbol{y}}=\boldsymbol{\theta}
+$$
+
+where the matrices $\mathbf{M}$ , $\mathbf{C}$ and $\mathbf{K}$ are the $T$ -periodic tangent system matrices that are employed in the modal analysis described in the next section.
+
+# 3. METHODS
+
+In this section, the four methods for modal analysis of structures with rotors are presented.
+
+# 3.1. Coleman approach
+
+The Coleman transformation requires identical degrees of freedom on each blade, and therefore, the equations of motion (equation (4)) in global coordinates were first transformed into substructure coordinates $y_{\mathrm{T}}$ . The transformation is
+
+$$
+\begin{array}{r l}&{\boldsymbol{y}=\mathrm{\mathbf{T}}\boldsymbol{y}_{\mathrm{T}}}\\ &{\mathbf{T}=\mathbf{diag}(\mathbf{I}_{N_{s}},\mathbf{T}_{\mathrm{r}},\mathbf{T}_{\mathrm{b1}},\mathbf{T}_{\mathrm{b2}},\mathbf{T}_{\mathrm{b3}})}\end{array}
+$$
+
+where $\mathbf{T}$ is a block diagonal time-variant matrix composed of the identity matrix $\mathbf{I}_{N_{\mathrm{s}}}$ sized by the number of degrees of freedom of the tower, the nacelle and the drivetrain, $\mathbf{T_{r}}$ transforms the degrees of freedom on the shaft and the hub into a hub centre frame and $\mathrm{T}_{\mathfrak{b}j}$ transforms the degrees of freedom on blade number $j=1,2,3$ into a local frame for blade $j$ . The triads were obtained in the periodic steady state, and thus, $\mathbf{T}$ is $T$ -periodic.
+
+The time-variant transformation into inertial frame coordinates $z$ is
+
+$$
+\begin{array}{r l}&{{\mathbf y}_{\mathrm{T}}={\mathbf B}\,z}\\ &{{\mathbf B}={\textbf d i a g}({\mathbf I}_{N_{\mathrm{s}}},{\mathbf B}_{\mathrm{r}},{\mathbf B}_{\mathrm{b}})}\end{array}
+$$
+
+where $\mathbf{B}_{\mathrm{r}}$ is a simple rotational transformation of the shaft and the hub and $\mathbf{B}_{\mathrm{b}}$ is the Coleman transformation introducing multiblade coordinates for a three-bladed rotor11,19 as
+
+$$
+\mathbf{B}_{\mathrm{b}}=\left[\begin{array}{l l l}{\mathbf{I}_{N_{\mathrm{b}}}}&{\mathbf{I}_{N_{\mathrm{b}}}\cos\psi_{1}}&{\mathbf{I}_{N_{\mathrm{b}}}\sin\psi_{1}}\\ {\mathbf{I}_{N_{\mathrm{b}}}}&{\mathbf{I}_{N_{\mathrm{b}}}\cos\psi_{2}}&{\mathbf{I}_{N_{\mathrm{b}}}\sin\psi_{2}}\\ {\mathbf{I}_{N_{\mathrm{b}}}}&{\mathbf{I}_{N_{\mathrm{b}}}\cos\psi_{3}}&{\mathbf{I}_{N_{\mathrm{b}}}\sin\psi_{3}}\end{array}\right]
+$$
+
+where $\psi_{j}=\varOmega t+2\pi(j-1)/3$ is the mean azimuth angle to blade number $j$ and $N_{\mathrm{b}}$ is the number of degrees of freedom on each blade. The inertial frame coordinate vector
+
+$$
+\boldsymbol{z}=\{y_{\mathrm{s}}^{\mathrm{T}}\,z_{\mathrm{r}}^{\mathrm{T}}\,a_{0}^{\mathrm{T}}\,a_{1}^{\mathrm{T}}\,b_{1}^{\mathrm{T}}\}^{\mathrm{T}}
+$$
+
+contains the untransformed coordinates for tower, nacelle and drivetrain ${\mathfrak{y}}_{\mathrm{s}}$ , the coordinates for shaft and hub $z_{\mathrm{r}}$ measured in a non-rotating frame aligned with the hub and the multiblade symmetric coordinates $\pmb{a}_{0}$ , cosine coordinates $\pmb{a}_{1}$ and sine coordinates $\pmb{b}_{1}$ . The details on how multiblade coordinates describe the motion of a wind turbine rotor in the inertial frame are discussed by Hansen.20,21
+
+The Coleman transformed equations were obtained by first inserting equation (5) into equation (4), then converting it to first order form and lastly introducing the inertial frame transformation in equation (6) as ${\displaystyle y_{\mathrm{T}2}=\mathrm{diag}({\bf B},{\bf B})z_{2}}$ where ${y}_{\mathrm{T}2}=\{{y}^{\mathrm{T}}\,{\dot{{y}}}^{\mathrm{T}}\}^{\mathrm{T}}$ and $z_{2}=\dot{\{z^{\mathrm{T}}\,\tilde{z}^{\mathrm{T}}\}}^{\mathrm{T}}$ are the state v ectors in substructure and ine rtial frames, respectively, with $\tilde{z}=\dot{z}+\bar{\omega}z$ and the constant matrix $\bar{\boldsymbol{\omega}}=\mathbf{B}^{-1}\mathbf{B}$ . The result is
+
+$$
+\begin{array}{r l}&{\dot{z}_{2}=\mathbf{A}_{\mathrm{B}}z_{2}}\\ &{\mathbf{A}_{\mathrm{B}}=\left[\mathbf{-}\mathbf{\bar{\omega}}-\mathbf{\bar{\omega}}\mathbf{\bar{\omega}}_{\mathrm{KB}}\quad\mathbf{-M}_{\mathrm{B}}^{-1}\mathbf{C}_{\mathrm{B}}-\mathbf{\bar{\omega}}\right]}\end{array}
+$$
+
+where $\mathbf{A}_{\mathrm{B}}$ is the Coleman transformed system matrix and
+
+$$
+\begin{array}{r l}&{\mathbf{M}_{\mathrm{B}}=\mathbf{B}^{-1}\mathbf{T}^{\mathrm{T}}\mathbf{M}\mathbf{T}\,\mathbf{B}}\\ &{\mathbf{C}_{\mathrm{B}}=\mathbf{B}^{-1}\mathbf{T}^{\mathrm{T}}(\mathbf{C}\,\mathbf{T}+2\,\mathbf{M}\,\dot{\mathbf{T}})\mathbf{B}}\\ &{\mathbf{K}_{\mathrm{B}}=\mathbf{B}^{-1}\mathbf{T}^{\mathrm{T}}(\mathbf{K}\,\mathbf{T}+\mathbf{C}\,\dot{\mathbf{T}}+\mathbf{M}\,\ddot{\mathbf{T}})\mathbf{B}}\end{array}
+$$
+
+are the Coleman transformed mass, damping/gyroscopic and stiffness matrices, respectively. If the system is isotropic, then $\mathbf{A}_{\mathrm{B}}$ is time-invariant, and a transient solution of equation (9) is
+
+$$
+z_{2}=\mathrm{e}^{\mathbf{A}_{\mathrm{B}}t}z_{2}(0)=\mathbf{V}\mathrm{e}^{\mathbf{A}t}q(0)
+$$
+
+where $\mathbf{A}$ is a diagonal matrix containing the eigenvalues of $\mathbf{A}_{\mathrm{B}}$ , $\mathbf{V}$ contains the corresponding eigenvectors as columns and $\pmb q(0)=\mathbf V^{-1}z_{2}(0)$ are the initial conditions in modal c oordinates. It is assumed th at all eigenvect ors are linearly independen t .
+
+The blade motion given in the inertial frame in equation (11) can be transformed back into the rotating frame using equation (6) as21
+
+$$
+\Gamma,i k=\mathrm{e}^{\sigma_{k}t}\left(A_{0,i k}\cos(\omega_{k}t+\varphi_{0,i k})+A_{\mathrm{BW},i k}\cos\left((\omega_{k}+\Omega)t+\varphi_{j}+\varphi_{\mathrm{BW},i k}\right)+A_{\mathrm{FW},i k}\cos\left((\omega_{k}-\Omega)t-\varphi_{j}+\varphi_{\mathrm{GW},i k}\right)\right),
+$$
+
+where $\varphi_{j}=2\pi(j-1)/3$ and $\sigma_{k}$ and $\omega_{k}$ pare the modal damping and frequency of mode number $k$ , respectively, given by the eigenvalue $\lambda_{k}=\sigma_{k}+\mathrm{i}\omega_{k}$ with $\mathrm{i}=\sqrt{-1}$ . The amplitudes for degree of freedom number $i$ were determined from the components of the eigenvector $\nu_{k}$ gi ven in multiblade coordinates of equation (8) as $A_{0,i k}=|a_{0,i k}|$ and
+
+$$
+\begin{array}{r l}&{A_{\mathrm{BW},i k}=\frac{1}{2}\big((\mathrm{Re}\,(a_{1,i k})+\mathrm{Im}\,(b_{1,i k}))^{2}+(\mathrm{Re}\,(b_{1,i k})-\mathrm{Im}\,(a_{1,i k}))^{2}\big)^{1/2}}\\ &{A_{\mathrm{FW},i k}=\frac{1}{2}\big((\mathrm{Re}\,(a_{1,i k})-\mathrm{Im}\,(b_{1,i k}))^{2}+(\mathrm{Re}\,(b_{1,i k})+\mathrm{Im}\,(a_{1,i k}))^{2}\big)^{1/2}}\end{array}
+$$
+
+where the subscripts 0, BW and FW denote symmetric, backward whirling and forward whirling motion, respectively.
+
+# 3.2. Classical Floquet analysis
+
+Floquet analysis enables the solution of the periodic equations of motion directly without an explicit transformation. Equation (4) is written in first order form
+
+$$
+\begin{array}{r l}&{\dot{\boldsymbol{y}}_{2}=\mathbf{A}\boldsymbol{y}_{2}}\\ &{\mathbf{A}=\left[\mathbf{-M}^{-1}\mathbf{K}\quad\mathbf{-M}^{-1}\mathbf{C}\right]}\end{array}
+$$
+
+where ${\mathbf y}_{2}=\{{\mathbf y}^{{\mathrm T}}\,\dot{{\mathbf y}}^{{\mathrm T}}\}^{{\mathrm T}}$ is the state vector and $\mathbf{A}$ is the $T$ -periodic system matrix.
+
+Floquet theory22 states that the solution to equation (15) is of the form
+
+$$
+\mathbf{\boldsymbol{y}}_{2}=\mathbf{\boldsymbol{U}}\mathbf{\boldsymbol{e}}^{\mathbf{\boldsymbol{\Lambda}}t}\mathbf{\boldsymbol{U}}^{-1}(0)\mathbf{\boldsymbol{y}}_{2}(0)
+$$
+
+where $\mathbf{U}$ is a $T$ -periodic matrix and $\mathbf{A}$ is a diagonal matrix. One way to construct this solution is to form a fundamental solution to equation (15) as
+
+$$
+\displaystyle\varphi=\bigl[\varphi_{1}\quad\varphi_{2}\quad.\ .\quad\varphi_{N}\bigr]
+$$
+
+over one period, $t\ \in\ [0;T]$ , where $N$ is the number of state variables, such that $\dot{\varphi}\;=\;{\bf A}\varphi$ . The monodromy matrix defined as
+
+$$
+\mathbf{C}=\boldsymbol{\varphi}^{-1}(0)\boldsymbol{\varphi}(T)
+$$
+
+contains all modal properties, which can be extracted from the eigenvalue decomposition
+
+$$
+\mathbf{C}=\mathbf{V}\mathbf{J}\mathbf{V}^{-1}
+$$
+
+where $\mathbf{V}$ contains the column eigenvectors $\nu_{k}$ of $\mathbf{C}$ , which are all assumed to be linearly independent and $\mathbf{J}$ is a diagonal matrix containing the eigenvalues $\rho_{k}$ of $\mathbf{C}$ , called the characteristic multipliers. The characteristic exponents $\lambda_{k}=\sigma_{k}+\mathrm{i}\omega_{k}$ contain the frequency $\omega_{k}$ and damping $\sigma_{k}$ and are related to the characteristic multipliers as $\rho_{k}=\exp(\lambda_{k}T)$ . Because the complex logarithm is not unique, the frequency is not determined uniquely, and the principal frequency $\omega_{\mathrm{p},k}$ and the damping $\sigma_{k}$ are defined from the characteristic multipliers as
+
+$$
+\begin{array}{c}{\displaystyle\sigma_{k}=\frac{1}{T}\ln(\vert\rho_{k}\vert)}\\ {\displaystyle\omega_{\mathrm{p},k}=\frac{1}{T}\arg(\rho_{k})}\end{array}
+$$
+
+where $\arg(\rho_{k})\in]-\pi;\pi]$ is implied, resulting in $\omega_{\mathrm{p},k}\in\big[-\Omega/2;\Omega/2\big]$ . Any integer multiple of the rotor speed can be added to the principal frequency to obtain a more physically meaningful frequency23,24
+
+$$
+\omega_{k}=\omega_{\mathrm{p},k}+j_{k}\Omega
+$$
+
+a choice that also affects the periodic modal matrix $\mathbf{U}$ in equation (16). This matrix $\mathbf{U}$ contains the periodic mode shapes uk and is given as24
+
+$$
+\pmb{u}_{k}=\varphi\nu_{k}\mathrm{e}^{-\lambda_{k}t}
+$$
+
+where the real part of $\lambda_{k}$ is given by equation (20) and the imaginary part of $\lambda_{k}$ is defined by equation (21) by selecting $j_{k}$ such that $\pmb{u}_{k}$ is as constant as possible for degrees of freedXom measured in the inertial frame.
+
+Introducing the Fourier transform of the periodic mode shape
+
+$$
+{\pmb u}_{k}=\sum_{j=-\infty}^{\infty}u_{j k}\mathrm{e}^{\mathrm{i}j\Omega t}
+$$
+
+the transient solution in equation (16) can be written as a sum of harmonic components
+
+$$
+y_{2}=\sum_{k=1}^{N}\sum_{j=-\infty}^{\infty}\mathcal{U}_{j k}\mathrm{e}^{(\sigma_{k}+\mathrm{i}(\omega_{k}+j\Omega))t}q_{k}(0)
+$$
+
+where $\begin{array}{r}{\pmb q(0)=\mathbf{U}^{-1}(0)\mathbf{y}_{2}(0).}\end{array}$ . Note that equation (12) is a special case of this expression for $j=-1,0,1$
+
+# 3.3. Implicit Floquet analysis
+
+The implicit Floquet method is here described based on the detailed description in Bauchau and Nikishkov,12 which focuses on computation of the characteristic multipliers from the state transition matrix $\Phi(T,0)$ . It can be defined in classical Floquet theory as
+
+$$
+\boldsymbol{\varphi}(T)=\boldsymbol{\Phi}(T,0)\,\boldsymbol{\varphi}(0)
+$$
+
+Using equation (18), the relationship between the state transition and monodromy matrices is derived as
+
+$$
+\Phi(T,0)=\varphi(0){\bf C}\,\varphi^{-1}(0)
+$$
+
+showing that $\Phi(T,0)$ and C have identical eigenvalues (characteristic multipliers), and their eigenvectors are related as $\pmb{\nu}_{k}=\pmb{\varphi}^{-1}(0)\pmb{w}_{k}$ , where $w_{k}$ represents the eigenvectors of $\Phi(T,0)$ .
+
+The key feature of the state transition matrix is that it defines the solution $y_{2}(T)=\Phi(T,0)y_{2}(0)$ for a time integration of the system equations (equation (15)) over one period $T$ with initial conditions ${\mathfrak{y}}_{2}(0)$ . Hence, without knowing the state transition matrix, it is possible to obtain the product of it with an arbitrary vector (the initial state vector) by integration of
+
+equation (15) over one period. The Arnoldi algorithm25 is a method to approximate the eigenvalues and the eigenvectors of a matrix, say $\Phi(T,0)$ , using only the matrix multiplication with $\Phi(T,0)$ to construct an $m$ -sized subspace
+
+$$
+\mathbf{P}=\left[p_{1}\quad p_{2}\quad\ldots\quad p_{m}\right]
+$$
+
+that satisfies the orthonormality condition
+
+$$
+\mathbf{P}^{\mathrm{T}}\mathbf{P}=\mathbf{I}
+$$
+
+and where the eigenvalues $\tilde{\rho}_{k}$ of the subspace projected state transition matrix
+
+$$
+\mathbf{H}=\mathbf{P}^{\mathrm{T}}\Phi(T,0)\mathbf{P}
+$$
+
+converge towards the eigenvalues $\rho_{k}$ of $\Phi(T,0)$ with the largest modulus as the size $m$ of the subspace increases. The subspace eigenvectors $\tilde{w}_{k}$ of $\mathbf{H}$ projected back to the full state space converge towards the eigenvectors $w_{k}$ of $\Phi(T,0)$ , i.e. $w_{k}\approx\mathbf{P}\tilde{w}_{k}$ . The Arnoldi algorithm proceeds as follows:
+
+Choose an arbitrary vector $\pmb{p}_{1}$ with $|p_{1}|=1$ for $n=1,2,\ldots,m$ $\pmb{a}:=\Phi(T,0)p_{n}$ (integration of equation (15) over $t\in[0;T])$ $\begin{array}{l}{b:=a}\\ {\mathrm{for~}j=1,2,\dotsc,n}\\ {\quad h_{j,n}:=p_{j}^{\operatorname{T}}a}\\ {\quad b:=b-h_{j,n}p_{j}}\end{array}$ end if $n