Merge pull request OpenFAST#3147 from JustinPorter88/modal-damping

Modal damping

Author: deslaughter

docs/source/user/api_change.rst

@@ -64,6 +64,9 @@ ElastoDyn                                     79       PBrIner(3)           200
 ElastoDyn                                     80       BlPIner(1)           28578         BlPIner(1)  - Pitch inertia of an undeflected blade, blade 1 (kg m^2)
 ElastoDyn                                     81       BlPIner(2)           28578         BlPIner(2)  - Pitch inertia of an undeflected blade, blade 2 (kg m^2)
 ElastoDyn                                     82       BlPIner(3)           28578         BlPIner(3)  - Pitch inertia of an undeflected blade, blade 3 (kg m^2) [unused for 2 blades]
+BeamDyn                                       10                            ------ Modal Damping [used only if damp_type=2] --------------------------------
+BeamDyn                                       11       n_modes              3             n_modes     - Number of modal damping coefficients (-)
+BeamDyn                                       12       zeta                 0.1, 0.2, 0.3 zeta        - Damping coefficients for mode 1 through n_modes
 ServoDyn                                      9        PitNeut(1)           0             PitNeut(1)  - Blade 1 neutral pitch position--pitch spring moment is zero at this position *[unused when* **PCMode>0** and **t>=TPCOn** *]*
 ServoDyn                                      10       PitNeut(2)           0             PitNeut(2)  - Blade 2 neutral pitch position--pitch spring moment is zero at this position *[unused when* **PCMode>0** and **t>=TPCOn** *]*
 ServoDyn                                      11       PitNeut(3)           0             PitNeut(3)  - Blade 3 neutral pitch position--pitch spring moment is zero at this position *[unused when* **PCMode>0** and **t>=TPCOn** *]* *[unused for 2 blades]*

docs/source/user/beamdyn/appendix.rst

@@ -16,6 +16,12 @@ OpenFAST+BeamDyn and stand-alone BeamDyn (static and dynamic) simulations all re
 :download:`(NREL 5MW static example) <examples/bd_primary_nrel_5mw.inp>`: This file includes information on the numerical-solution parameters (e.g., numerical damping, quadrature rules), and the geometric definition of the beam reference line via "members" and "key points".  This file also specifies the "blade input file."
 
 2) BeamDyn blade input file :download:`(NREL 5MW example) <examples/nrel_5mw_blade.inp>`: 
+   This file specifies the blade sectional properties at various stations along the blade.
+   The file includes stiffness and mass matrices at each station, as well as damping parameters.
+   Note that the example file uses stiffness-proportional damping (damp_flag = 1). For modal
+   damping (damp_flag = 2), the n_modes parameter should be set to a non-zero value and 
+   followed by the corresponding modal damping ratios (zeta values) represented as the
+   fraction of critical damping.
 
 Stand-alone BeamDyn simulation also require a driver input file; we list here examples for static and dynamic simulations:
 
@@ -56,3 +62,21 @@ outputs are expressed in one of the following three coordinate systems:
    :align: center
 
    BeamDyn Output Channel List
+
+.. note::
+
+   **New Output Channels (v5.0 and later):**
+   
+   BeamDyn now includes additional output channels for applied loads mapped to the root node.
+   These channels provide the total applied loads (both distributed and point loads) resolved
+   at the root of the blade, expressed in both the root coordinate system (r) and global 
+   inertial frame (g):
+   
+   - **RootAppliedFxr, RootAppliedFyr, RootAppliedFzr**: Applied force components in r-frame
+   - **RootAppliedMxr, RootAppliedMyr, RootAppliedMzr**: Applied moment components in r-frame
+   - **RootAppliedFxg, RootAppliedFyg, RootAppliedFzg**: Applied force components in g-frame
+   - **RootAppliedMxg, RootAppliedMyg, RootAppliedMzg**: Applied moment components in g-frame
+   
+   These outputs are useful for understanding the total aerodynamic and other external loads
+   acting on the blade, particularly when diagnosing load imbalances or validating force 
+   distributions.

docs/source/user/beamdyn/examples/nrel_5mw_blade.inp

@@ -1,13 +1,16 @@
  ------- BEAMDYN V1.00.* INDIVIDUAL BLADE INPUT FILE --------------------------
  NREL 5MW Blade
- ---------------------- BLADE PARAMETERS --------------------------------------
-49   station_total    - Number of blade input stations (-)
- 1   damp_flag        - Damping flag: 0: no damping; 1: damped
-  ---------------------- DAMPING COEFFICIENT------------------------------------
+------ Blade Parameters --------------------------------------------------------
+49                      station_total    - Number of blade input stations (-)
+1                       damp_type        - Damping type (switch) {0: none, 1: stiffness-proportional, 2: modal}
+------ Stiffness-Proportional Damping [used only if damp_type=1] ---------------
    mu1        mu2        mu3        mu4        mu5        mu6
    (-)        (-)        (-)        (-)        (-)        (-)
-1.0E-03    1.0E-03    1.0E-03    0.0014      0.0022     0.0022
- ---------------------- DISTRIBUTED PROPERTIES---------------------------------
+1.0E-03    1.0E-03    1.0E-03     0.0014      0.0022     0.0022
+------ Modal Damping [used only if damp_type=2] --------------------------------
+4                        n_modes    - Number of modal damping coefficients (-)
+0.01 0.02 0.03 0.04      zeta       - Damping coefficients for mode 1 through n_modes
+------ Distributed Properties --------------------------------------------------
   0.000000
    9.729480E+08    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00
    0.000000E+00    9.729480E+08    0.000000E+00    0.000000E+00    0.000000E+00    0.000000E+00

docs/source/user/beamdyn/input_files.rst

@@ -215,15 +215,22 @@ contents of the BeamDyn input file (useful for debugging errors in the
 input file).
 
 The ``QuasiStaticInit`` flag indicates if BeamDyn should perform a quasi-static
-solution at initialization to better initialize its states. In general, this should 
-be set to true for better numerical performance (it reduces startup transients).
+solution at initialization to better initialize its states. This option is only 
+available when coupled to OpenFAST with dynamic analysis enabled. When set to ``TRUE``,
+BeamDyn will perform a steady-state startup (SSS) solve that includes centripetal
+accelerations to reduce startup transients and improve numerical performance. This
+is particularly useful for rotating blade simulations where the initial conditions
+would otherwise cause spurious transients. The keyword "DEFAULT" sets this to ``FALSE``.
 
 ``rhoinf`` specifies the numerical damping parameter (spectral radius
 of the amplification matrix) in the range of :math:`[0.0,1.0]` used in
 the generalized-\ :math:`\alpha` time integrator implemented in BeamDyn
-for dynamic analysis. For ``rhoinf = 1.0``, no
-numerical damping is introduced and the generalized-\ :math:`\alpha`
-scheme is identical to the Newmark scheme; for
+for dynamic analysis. **Note: This parameter is only used when BeamDyn is run
+in loose coupling mode (e.g., stand-alone with the driver or when loose coupling 
+is explicitly selected in OpenFAST). When tight coupling is used in OpenFAST, 
+time integration is handled by the glue code and this parameter is ignored.**
+For ``rhoinf = 1.0``, no numerical damping is introduced and the 
+generalized-\ :math:`\alpha` scheme is identical to the Newmark scheme; for
 ``rhoinf = 0.0``, maximum numerical damping is
 introduced. Numerical damping may help produce numerically stable
 solutions.
@@ -245,9 +252,11 @@ between two input stations into “Refine factor” of segments. The keyword
 This entry is not used in Gauss quadrature.
 
 ``N_Fact`` specifies a parameter used in the modified Newton-Raphson
-scheme. If ``N_Fact = 1`` a full Newton
-iteration scheme is used, i.e., the global tangent stiffness matrix is
-computed and factorized at each iteration; if
+scheme. **Note: This parameter is only used when BeamDyn is run in loose coupling 
+mode (stand-alone with the driver or when loose coupling is explicitly selected in 
+OpenFAST). When tight coupling is used in OpenFAST, this parameter is ignored.**
+If ``N_Fact = 1`` a full Newton iteration scheme is used, i.e., the global tangent 
+stiffness matrix is computed and factorized at each iteration; if
 ``N_Fact > 1`` a modified Newton iteration
 scheme is used, i.e., the global stiffness matrix is computed and
 factorized every ``N_Fact`` iterations within each time step. The
@@ -265,15 +274,19 @@ load steps as opposed to one single large load step which may cause divergence o
 Newton-Raphson scheme. The keyword “DEFAULT” sets ``load_retries = 20``.
 
 ``NRMax`` specifies the maximum number of iterations per time step in
-the Newton-Raphson scheme. If convergence is not reached within this
-number of iterations, BeamDyn returns an error message and terminates
-the simulation. The keyword “DEFAULT” sets
+the Newton-Raphson scheme. **Note: This parameter is only used when BeamDyn is run 
+in loose coupling mode (stand-alone with the driver or when loose coupling is explicitly 
+selected in OpenFAST). When tight coupling is used in OpenFAST, this parameter is ignored.**
+If convergence is not reached within this number of iterations, BeamDyn returns an error 
+message and terminates the simulation. The keyword "DEFAULT" sets
 ``NRMax = 10``.
 
 ``Stop_Tol`` specifies a tolerance parameter used in convergence
 criteria of a nonlinear solution that is used for the termination of the
-iteration. The keyword “DEFAULT” sets
-``Stop_Tol = 1.0E-05``. Please refer to
+iteration. **Note: This parameter is only used when BeamDyn is run in loose coupling 
+mode (stand-alone with the driver or when loose coupling is explicitly selected in 
+OpenFAST). When tight coupling is used in OpenFAST, this parameter is ignored.**
+The keyword "DEFAULT" sets ``Stop_Tol = 1.0E-05``. Please refer to
 :numref:`convergence-criterion` for more details.
 
 ``tngt_stf_fd`` is a boolean that sets the flag to compute the tangent stiffness
@@ -468,18 +481,26 @@ Blade Parameters
 ``Station_Total`` specifies the number cross-sectional stations along
 the blade axis used in the analysis.
 
-``Damp_Type`` specifies if structural damping is considered in the
-analysis. If ``Damp_Type = 0``, then no damping is
-considered in the analysis and the six damping coefficient in the next
-section will be ignored. If ``Damp_Type = 1``, structural
-damping will be included in the analysis.
+``damp_type`` specifies the type of structural damping to be used in the
+analysis. There are three options:
 
-Damping Coefficient
-~~~~~~~~~~~~~~~~~~~
+- ``damp_type = 0``: No damping is considered in the analysis. The damping 
+  coefficients in the following sections will be ignored.
+- ``damp_type = 1``: Stiffness-proportional damping is applied. 
+  The six damping coefficients (``beta``) in the Stiffness-Proportional Damping 
+  section are used to scale the 6x6 stiffness matrix at each cross section.
+- ``damp_type = 2``: Modal damping is applied. The modal fractions of critical damping
+  (``zeta``) for the first ``n_modes`` modes are used. BeamDyn internally computes 
+  the modal properties and applies damping in the modal coordinates.
+
+Stiffness-Proportional Damping Coefficients
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-This section specifies six damping coefficients, :math:`\mu_{ii}` with
+This section specifies six damping coefficients, :math:`\mu_{ii}` (also called ``beta``) with
 :math:`i \in [1,6]`, for six DOFs (three translations and three
-rotations). Viscous damping is implemented in BeamDyn where the damping
+rotations). These coefficients are only used when ``damp_flag = 1``.
+
+Viscous damping is implemented in BeamDyn where the damping
 forces are proportional to the strain rate. These are
 stiffness-proportional damping coefficients, whereby the
 :math:`6\times6` damping matrix at each cross section is scaled from the
@@ -510,6 +531,53 @@ coefficient matrix defined as
        0 & 0 & 0 & 0 & 0 & \mu_{66} \\
    \end{bmatrix}
 
+Modal Damping Coefficients
+~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This section specifies modal damping parameters and is only used when ``damp_flag = 2``.
+
+``n_modes`` specifies the number of modes for which modal damping coefficients
+are provided. BeamDyn will compute the natural frequencies and mode shapes of the 
+blade and apply damping in the modal coordinates.
+
+``zeta`` is an array of ``n_modes`` modal damping ratios, one for each mode. 
+Each value should typically be between 0.0 (no damping) and 1.0 (critical damping),
+though higher values are permitted. Common values for composite blade structures
+are less than 0.05 (less than 5% of critical damping). The damping
+ratios are applied to modes 1 through ``n_modes`` in order of increasing frequency.
+
+After ``n_modes`` damping values are assigned that grow proportional to the modal
+natural frequency. This proportionality is assigned to match the ``zeta`` value
+of last mode with prescribed damping.
+
+If modal damping is selected, BeamDyn calculates nodal damping forces based on the node velocities, 
+rotated to the initial node orientation, and the mode shape after quasi-static initialization has been performed,
+if it was requested. These nodal damping forces are then transformed back to the current node orientation.
+
+When using modal damping, ``PitchDOF`` needs to be set to ``TRUE`` in ElastoDyn to give BeamDyn
+a correct root pitch velocity. Otherwise, blade pitching will give spurious damping behavior.
+
+Recommendations:
+
+- It is recommended to stop inputting zeta values before reaching the first axial mode (typically around 18 modes).
+  Users may experiment with including more or fewer modes to observe the effect on results.
+- Avoid prescribing a final zeta value of 1.0 (e.g., do not specify 18 modes with realistic zeta values followed by a 19th mode with zeta=1.0), as this can significantly degrade result quality.
+- When attempting to match stiffness-proportional damping (:math:`\mu`), the OpenFAST toolbox may fail to
+  provide reliable damping values matched with mode numbers once some modes become critically damped.
+  Reducing the number of modes (e.g., from 40 to 30) can help if higher modes are indexed incorrectly.
+- In some cases, axial loads appear to be driven by the axial motion of non-axial modes.
+- Verify that the modal frequencies and modal participation factors output in ``InputFile.BD.R*.B*.modes.csv`` file are consistent with the expected modes.
+  This file is only output when modal damping is used.
+  The number of output modes in this file is equal to the degrees of freedom of the blade.
+  However, modes tend to become numerical artifacts past the first few modes of a given type.
+
+Notes:
+
+- Modal damping is experimental in OpenFAST 5.0, so should be used with care.
+- Modal damping is mainly supported for tight coupling with OpenFAST.
+- In some cases, when used with loose coupling or BeamDyn Driver, a much smaller time step may be needed to get reaction forces to converge with the tight coupling case and/or stiffness-proportional damping (:math:`\mu`) case. 
+- Output forces internal forces from BeamDyn at quadrature points do not capture contributions from the modal damping forces.
+
 Distributed Properties
 ~~~~~~~~~~~~~~~~~~~~~~
 

docs/source/user/beamdyn/introduction.rst

@@ -5,14 +5,14 @@ Introduction
 
 BeamDyn is a time-domain structural-dynamics module for slender
 structures created by the National Renewable Energy Laboratory (NREL)
-through support from the U.S. Department of Energy Wind and Water Power
+through support from the U.S. Department of Energy Wind and Water Power
 Program and the NREL Laboratory Directed Research and Development (LDRD)
 program through the grant “High-Fidelity Computational Modeling of
-Wind-Turbine Structural Dynamics”, see References :cite:`Wang:SFE2013,Wang:GEBT2013,Wang:GEBT2014,Wang:2015`.
+Wind-Turbine Structural Dynamics”, see References :cite:`Wang:SFE2013,Wang:GEBT2013,Wang:GEBT2014,Wang:2015`.
 The module has been coupled into the FAST aero-hydro-servo-elastic wind
 turbine multi-physics engineering tool where it used to model blade
 structural dynamics. The BeamDyn module follows the requirements of the
-FAST modularization framework, see References :cite:`Jonkman:2013`;
+FAST modularization framework, see References :cite:`Jonkman:2013`;
 :cite:`Sprague:2013,Sprague:2014,website:FASTModularizationFramework`,
 couples to FAST version 8, and provides new capabilities for modeling
 initially curved and twisted composite wind turbine blades undergoing
@@ -64,7 +64,7 @@ outputs the blade displacements, velocities, and accelerations along the
 beam length, which are used by AeroDyn to calculate the local
 aerodynamic loads (distributed along the length) that are used as inputs
 for BeamDyn. In addition, BeamDyn can calculate member internal reaction
-loads, as requested by the user. Please refers to Figure [fig:FlowChart]
+loads, as requested by the user. Please refers to Figure [fig:FlowChart]
 for the coupled interactions between BeamDyn and other modules in FAST.
 When coupled to FAST, BeamDyn replaces the more simplified blade
 structural model of ElastoDyn that is still available as an option, but
@@ -100,12 +100,12 @@ loads specified at FE nodes, or a combination of the two. When BeamDyn
 is coupled to FAST, the blade analysis node discretization may be
 independent between BeamDyn and AeroDyn.
 
-This document is organized as follows. Section :ref:`running-beamdyn` details how to
+This document is organized as follows. Section :ref:`running-beamdyn` details how to
 obtain the BeamDyn and FAST software archives and run either the
 stand-alone version of BeamDyn or BeamDyn coupled to FAST.
-Section :ref:`bd-input-files` describes the BeamDyn input files.
-Section :ref:`bd-output-files` discusses the output files generated by
-BeamDyn. Section :ref:`beamdyn-theory` summarizes the BeamDyn theory.
-Section :ref:`bd-future-work` outlines potential future work. Example input
-files are shown in Appendix :numref:`bd_input_files`.
-A summary of available output channels is found in Appendix :ref:`app-output-channel`.
+Section :ref:`bd-input-files` describes the BeamDyn input files.
+Section :ref:`bd-output-files` discusses the output files generated by
+BeamDyn. Section :ref:`beamdyn-theory` summarizes the BeamDyn theory.
+Section :ref:`bd-future-work` outlines potential future work. Example input
+files are shown in Appendix :numref:`bd_input_files`.
+A summary of available output channels is found in Appendix :ref:`app-output-channel`.

docs/source/user/beamdyn/output_files.rst

@@ -26,7 +26,7 @@ Summary File
 In stand-alone mode, BeamDyn generates a summary file with the naming
 convention, ``InputFile.sum`` if the ``SumPrint`` parameter is set
 to TRUE. When coupled to FAST, the summary file is named
-``InputFile.BD.sum``. This file summarizes key information about the
+``InputFile.BD.R*.B*.sum.yaml``. This file summarizes key information about the
 simulation, including:
 
 -  Blade mass.