diff --git a/docs/source/user/api_change.rst b/docs/source/user/api_change.rst index 80a973a4ec..4f9269b09b 100644 --- a/docs/source/user/api_change.rst +++ b/docs/source/user/api_change.rst @@ -62,6 +62,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]* diff --git a/docs/source/user/beamdyn/appendix.rst b/docs/source/user/beamdyn/appendix.rst index 3e282ce5ae..3fabb468f7 100644 --- a/docs/source/user/beamdyn/appendix.rst +++ b/docs/source/user/beamdyn/appendix.rst @@ -16,6 +16,11 @@ OpenFAST+BeamDyn and stand-alone BeamDyn (static and dynamic) simulations all re :download:`(NREL 5MW static example) `: 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) `: + 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). Stand-alone BeamDyn simulation also require a driver input file; we list here examples for static and dynamic simulations: @@ -56,3 +61,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. diff --git a/docs/source/user/beamdyn/examples/nrel_5mw_blade.inp b/docs/source/user/beamdyn/examples/nrel_5mw_blade.inp index 24e2192f2a..0a5ebced03 100644 --- a/docs/source/user/beamdyn/examples/nrel_5mw_blade.inp +++ b/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 diff --git a/docs/source/user/beamdyn/input_files.rst b/docs/source/user/beamdyn/input_files.rst index 7f9741eea1..04b55ace47 100644 --- a/docs/source/user/beamdyn/input_files.rst +++ b/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 damping coefficients + (``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,32 @@ 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 in the range of 0.01 to 0.05 (1% to 5% of critical damping). The damping +ratios are applied to modes 1 through ``n_modes`` in order of increasing frequency. + +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. + +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. + Distributed Properties ~~~~~~~~~~~~~~~~~~~~~~ diff --git a/docs/source/user/beamdyn/introduction.rst b/docs/source/user/beamdyn/introduction.rst index 95e54cfc36..9cee56289b 100644 --- a/docs/source/user/beamdyn/introduction.rst +++ b/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`. diff --git a/docs/source/user/beamdyn/running_bd.rst b/docs/source/user/beamdyn/running_bd.rst index a2885315fc..2f0d767aba 100644 --- a/docs/source/user/beamdyn/running_bd.rst +++ b/docs/source/user/beamdyn/running_bd.rst @@ -14,41 +14,10 @@ There are two forms of the BeamDyn software to choose from: stand-alone and coupled to the FAST simulator. Although the user may not necessarily need both forms, he/she would likely need to be familiar with and run the stand-alone model if building a model of the blade from scratch. The -stand-alone version is also helpful for model troubleshooting, even if +stand-alone driver is also helpful for model troubleshooting, even if the goal is to conduct aero-hydro-servo-elastic simulations of onshore/offshore wind turbines within FAST. -Stand-Alone BeamDyn Archive -~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -Users can download the stand-alone BeamDyn archive from our Web server -at https://nwtc.nrel.gov/BeamDyn. The file has a name similar to -``BD_v1.00.00a.exe``, but may have a different version number. The user -can then download the self-extracting archive (*.exe*) to expand the -archive into a folder he/she specifies. - -The archive contains the ``bin``, ``CertTest``, ``Compiling``, -``Docs``, and ``Source`` folders. The ``bin`` folder includes the -main executable file, ``BeamDyn_Driver.exe``, which is used to execute -the stand-alone BeamDyn program. The ``CertTest`` folder contains a -collection of sample BeamDyn input files and driver input files that can -be used as templates for the user’s own models. This document may be -found in the ``Docs`` folder. The ``Compiling`` folder contains files -for compiling the stand-alone ``BeamDyn_v1.00.00.exe`` file with either -Visual Studio or gFortran. The Fortran source code is located in the -``Source`` folder. - -FAST Archive -~~~~~~~~~~~~ - -Download the FAST archive, which includes BeamDyn, from our Web server -at https://nwtc.nrel.gov/FAST8. The file has a name similar to -``FAST_v8.12.00.exe``, but may have a different version number. Run the -downloaded self-extracting archive (``.exe``) to expand the archive into a -user-specified folder. The FAST executable file is located in the -archive’s ``bin`` folder. An example model using the NREL 5-MW -reference turbine is located in the ``CertTest`` folder. - Running BeamDyn --------------- @@ -56,7 +25,7 @@ Running the Stand-Alone BeamDyn Program ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The stand-alone BeamDyn program, ``BeamDyn_Driver.exe``, simulates static -and dynamic responses of the user’s input model, without coupling to +and dynamic responses of the user's input model, without coupling to FAST. Unlike the coupled version, the stand-alone software requires the use of a driver file in addition to the primary and blade BeamDyn input files. This driver file specifies inputs normally provided to BeamDyn by @@ -78,12 +47,12 @@ file, as described in Section :ref:`driver-input-file`. Running BeamDyn Coupled to FAST ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -Run the coupled FAST software from a DOS command prompt by typing, for +Run the coupled FAST software from a command prompt by typing, for example: .. code-block:: bash - >FAST_Win32.exe Test26.fst + >openfast.exe Test26.fst where ``Test26.fst`` is the name of the primary FAST input file. This input file has a feature switch to enable or disable the BeamDyn diff --git a/docs/source/user/beamdyn/theory.rst b/docs/source/user/beamdyn/theory.rst index 78d14e6936..f6c590eaea 100644 --- a/docs/source/user/beamdyn/theory.rst +++ b/docs/source/user/beamdyn/theory.rst @@ -605,7 +605,32 @@ where :math:`{\underline{\underline{O}}}_{12}` and matrices of :math:`\mathcal{{\underline{\underline{O}}}}` and :math:`\mathcal{{\underline{\underline{G}}}}` as :math:`{\underline{\underline{C}}}_{12}` in Eq. :eq:`E1-PartC`. +Modal Damping +------------- +In addition to the stiffness-proportional viscous damping described above, BeamDyn +also supports modal damping. When modal damping is selected (``damp_flag = 2``), +BeamDyn computes the natural frequencies and mode shapes of the blade and applies +damping in the modal coordinates. + +The modal damping approach constructs a damping matrix :math:`{\underline{\underline{C}}}_{modal}` +based on user-specified modal damping ratios :math:`\zeta_i` for modes :math:`i = 1, 2, \ldots, n_{modes}`. +The modal damping matrix is defined in terms of the modal properties: + +.. math:: + :label: ModalDamping + + {\underline{\underline{C}}}_{modal} = \sum_{i=1}^{n_{modes}} 2 \zeta_i \omega_i \underline{\phi}_i \underline{\phi}_i^T + +where :math:`\omega_i` is the natural frequency of mode :math:`i`, :math:`\underline{\phi}_i` +is the corresponding mode shape (mass-normalized eigenvector), and :math:`\zeta_i` is the +modal damping ratio for mode :math:`i`. This damping matrix is then transformed to the +physical coordinates and applied during the time integration in the same manner as the +stiffness-proportional damping. + +The advantage of modal damping is that it allows for mode-specific damping levels, which +can be more physically representative of composite blade structures where different modes +may experience different levels of damping due to various energy dissipation mechanisms. .. _convergence-criterion: Convergence Criterion and Generalized-\ :math:`\alpha` Time Integrator @@ -637,7 +662,8 @@ Time integration is performed using the generalized-\ :math:`\alpha` scheme in BeamDyn, which is an unconditionally stable (for linear systems), second-order accurate algorithm. The scheme allows for users to choose integration parameters that introduce high-frequency numerical -dissipation. More details regarding the generalized-\ :math:`\alpha` +dissipation. The amount of numerical dissipation is controlled by the user-specified +spectral radius at infinity, :math:`\rho_{\infty}`. More details regarding the generalized-\ :math:`\alpha` method can be found in :cite:`Chung-Hulbert:1993,Bauchau:2010`. Calculation of Reaction Loads @@ -651,14 +677,14 @@ mode), the reaction loads at the root are needed to satisfy equality of the governing equations. The reaction loads at the root are also the loads passing from blade to hub in a full turbine analysis. -The governing equations in Eq. :eq:`GovernGEBT-1-2` can be recast in a compact form +The governing equations in Eq. :eq:`GovernGEBT-1-2` can be recast in a compact form .. math:: :label: CompactGovern {\underline{\mathcal{F}}}^I - {\underline{\mathcal{F}}}^{C\prime} + {\underline{\mathcal{F}}}^D = {\underline{\mathcal{F}}}^{ext} -with all the vectors defined in Section [sec:LinearProcess]. At the +with all the vectors defined in Section [sec:LinearProcess]. At the blade root, the governing equation is revised as .. math:: @@ -668,14 +694,14 @@ blade root, the governing equation is revised as where :math:`{\underline{\mathcal{F}}}^R = \left[ {\underline{F}}^R~~~{\underline{M}}^R\right]^T` is the reaction force vector and it can be solved from -Eq. :eq:`CompactGovernRoot` given that the motion fields are known at this +Eq. :eq:`CompactGovernRoot` given that the motion fields are known at this point. Calculation of Blade Loads -------------------------- BeamDyn can also calculate the blade loads at each finite element node -along the blade axis. The governing equation in Eq. :eq:`CompactGovern` are +along the blade axis. The governing equation in Eq. :eq:`CompactGovern` are recast as .. math:: diff --git a/modules/beamdyn/src/BeamDyn.f90 b/modules/beamdyn/src/BeamDyn.f90 index 6e245c482f..8c71621ec8 100644 --- a/modules/beamdyn/src/BeamDyn.f90 +++ b/modules/beamdyn/src/BeamDyn.f90 @@ -24,6 +24,8 @@ MODULE BeamDyn IMPLICIT NONE + INTEGER, PARAMETER :: FEKi = R8Ki ! Define the kind to be used for FEM/eigenanalysis + #ifndef UNIT_TEST PRIVATE @@ -188,7 +190,12 @@ SUBROUTINE BD_Init( InitInp, u, p, x, xd, z, OtherState, y, MiscVar, Interval, I QuasiStaticInitialized = .FALSE. ENDIF - + !................................. + ! Calculation of modal damping here + + IF(p%damp_flag .EQ. 2) THEN + call Init_ModalDamping(x, OtherState, p, MiscVar, ErrStat2, ErrMsg2); if (Failed()) return + ENDIF !................................. ! initialization of output mesh values (used for initial guess to AeroDyn) @@ -1102,7 +1109,16 @@ subroutine SetParameters(InitInp, InputFileData, p, OtherState, ErrStat, ErrMsg) ! Physical damping flag and 6 damping coefficients !............................................... p%damp_flag = InputFileData%InpBl%damp_flag - p%beta = InputFileData%InpBl%beta + select case (p%damp_flag) + case (0) ! No damping + case (1) ! Stiffness-proportional damping + p%beta = InputFileData%InpBl%beta + case (2) ! Modal damping + p%zeta = InputFileData%InpBl%zeta + case default + call SetErrStat(ErrID_Fatal, "Invalid value for physical damping flag in input file.", ErrStat, ErrMsg, RoutineName ) + return + end select !............................................... ! set parameters for File I/O data: @@ -1787,6 +1803,200 @@ END SUBROUTINE cleanup END SUBROUTINE Init_ContinuousStates +!----------------------------------------------------------------------------------------------------------------------------------- +!> This routine initializes modal damping. +SUBROUTINE Init_ModalDamping(x, OtherState, p, m, ErrStat, ErrMsg) + + TYPE(BD_ContinuousStateType), INTENT(IN ) :: x !< Continuous states at t on input at t + dt on output + type(BD_OtherStateType), INTENT(IN ) :: OtherState !< Global rotations are stored in otherstate + TYPE(BD_ParameterType), INTENT(INOUT) :: p !< Parameters, output modal damping matrix in original frame here + TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< misc/optimization variables + INTEGER(IntKi), INTENT( OUT) :: ErrStat + CHARACTER(*), INTENT( OUT) :: ErrMsg + + CHARACTER(*), PARAMETER :: RoutineName = 'Init_ModalDamping' + INTEGER(IntKi) :: ErrStat2 + CHARACTER(ErrMsgLen) :: ErrMsg2 + INTEGER(IntKi) :: nDOF + INTEGER(IntKi) :: j, k ! looping indexing variable + INTEGER(IntKi) :: numZeta ! number of damping values + REAL(R8Ki) :: Zj ! diagonal element of the modal damping matrix + REAL(R8Ki), ALLOCATABLE :: eigenvectors(:, :) ! mode shapes + REAL(R8Ki), ALLOCATABLE :: omega(:) ! modal frequencies (rad/s) + REAL(R8Ki), ALLOCATABLE :: phiT_M(:, :) ! mode shapes transpose times mass matrix + REAL(R8Ki), ALLOCATABLE :: phi0T_M_phi0(:, :) ! normalization calculation of mass matrix + REAL(R8Ki), ALLOCATABLE :: StifK(:,:) ! Copy of stiffness matrix for eigenanalysis (modified during solve) + REAL(R8Ki), ALLOCATABLE :: MassM(:,:) ! Copy of mass matrix for eigenanalysis (modified during solve) + real(R8Ki) :: NodeRot(3, 3) + + ErrStat = ErrID_None + ErrMsg = '' + + ! TODO : Take actual user input for zeta + ! zeta is fraction of critical damping. + ! zeta = (/ 0.1d0, 0.3d0, 0.15d0, 0.45d0, 0.2d0, 0.6d0, 0.7d0, 0.8d0, 0.9d0, 1.0d0 /) + ! zeta = (/ 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0 /) + + ! ! Zeta for Damped_Beam tests + ! zeta = (/ 0.001d0, 0.003d0, 0.0015d0, 0.0045d0, 0.002d0, 0.006d0, 0.007d0, 0.008d0, 0.009d0, 0.010d0 /) + + ! Zeta for IEA 22 test + ! zeta = (/ 0.00481431, 0.00501452, 0.01319441, 0.01365751, 0.02704056, 0.02928076, & + ! 0.03943607, 0.0133448, 0.05628846, 0.05935411, 0.01490789, 0.08727722, & + ! 0.08508602, 0.01754962, 0.12212305, 0.12607313, 0.02324255, 0.03571866, & + ! 0.13962916, 0.00848232, 0.16332203, 0.0405824, 0.28314572, 0.05114984, & + ! 0.40216838, 0.3507376, 0.022072, 0.0540725, 0.06678422, 0.03625773, & + ! 0.0484979, 0.8516855, 0.98869347, 0.05778878, 0.06610646, 0.0764022, & + ! 0.08499995, 0.09210838, 0.13037341, 0.13997288 /) + + numZeta = size(p%zeta) + + ! 0. Setup quadrature points + CALL BD_QuadraturePointData(p, x, m) + + ! 1. Generates K, M Matrices + ! These go into 'm%StifK' and 'm%MassM' + CALL BD_GenerateDynamicElementGA2( x, OtherState, p, m, .TRUE.) + + ! 2. Copy lines from 'BD_CalcForceAcc' for M, K -> 2D and apply Boundary conditions + ! Full mass matrix (n_dof, n_dof) + m%LP_MassM = reshape(m%MassM, [p%dof_total, p%dof_total]) + + ! Mass matrix for free nodes + m%LP_MassM_LU = m%LP_MassM(7:p%dof_total, 7:p%dof_total) + + ! Full stiffness matrix (n_dof, n_dof) + m%LP_StifK = reshape(m%StifK, [p%dof_total, p%dof_total]) + + ! Stiffness matrix for free nodes + m%LP_StifK_LU = m%LP_StifK(7:p%dof_total, 7:p%dof_total) + + ! 3. Do eigenanalysis + ! For now, calculate all eigenpairs + nDOF = p%dof_total - 6 + + ! Allocate eigenvector matrix and eigenvalue arrays + call AllocAry(eigenvectors, nDOF, nDOF, 'eigenvectors', ErrStat2, ErrMsg2); if (Failed()) return + call AllocAry(omega, nDOF, 'omega', ErrStat2, ErrMsg2); if (Failed()) return + call AllocAry(StifK, nDOF, nDOF, 'StifK', ErrStat2, ErrMsg2); if (Failed()) return + call AllocAry(MassM, nDOF, nDOF, 'MassM', ErrStat2, ErrMsg2); if (Failed()) return + + ! EigenSolve modifies the input matrices, so make copies before calling + StifK = m%LP_StifK_LU + MassM = m%LP_MassM_LU + call EigenSolve(StifK, MassM, nDOF, .TRUE., eigenvectors, omega, ErrStat2, ErrMsg2); if (Failed()) return + + ! Mass-normalize the mode shapes + call AllocAry(phi0T_M_phi0, nDOF, nDOF, 'phi0T_M_phi0', ErrStat2, ErrMsg2); if (Failed()) return + + phi0T_M_phi0 = matmul(transpose(eigenvectors), matmul(m%LP_MassM_LU, eigenvectors)) + + do j = 1, nDOF + eigenvectors(:, j) = eigenvectors(:, j) / sqrt(phi0T_M_phi0(j, j)) + end do + + ! 4. Generate damping matrix in original frame + call AllocAry(phiT_M, nDOF, nDOF, 'phiT_M', ErrStat2, ErrMsg2); if (Failed()) return + + phiT_M = matmul(transpose(eigenvectors), m%LP_MassM_LU) ! after normalization + + call AllocAry(p%ModalDampingMat, nDOF, nDOF, 'p%ModalDampingMat', ErrStat2, ErrMsg2); if (Failed()) return + + do j = 1, nDOF + + if( j <= numZeta) then + Zj = 2.0_R8Ki * omega(j) * p%zeta(j) + else + ! Stiffness proportional damping is used past the last prescribed value + ! at a rate equal to the last prescribed value. + Zj = 2.0_R8Ki * omega(j) * (p%zeta(numZeta) * omega(j) / omega(numZeta)) + endif + + p%ModalDampingMat(j, :) = Zj * phiT_M(j, :) + end do + + p%ModalDampingMat = matmul(transpose(phiT_M), p%ModalDampingMat) + + ! Apply the rotation of q here. When the actual dynamics are at the same position + ! as this, then this cancels with a rotation applied at the modal damping force. + do j = 2, p%node_total + + ! Loop over the nodes that apply to the damping matrix, so don't include the root node. + call BD_CrvMatrixR(x%q(4:6, j), NodeRot) + + k = (j-2)*6 + p%ModalDampingMat(:, k+1:k+3) = matmul(p%ModalDampingMat(:, k+1:k+3), NodeRot) + p%ModalDampingMat(:, k+4:k+6) = matmul(p%ModalDampingMat(:, k+4:k+6), NodeRot) + end do + do j = 2, p%node_total + + ! Loop over the nodes that apply to the damping matrix, so don't include the root node. + call BD_CrvMatrixR(x%q(4:6, j), NodeRot) + + k = (j-2)*6 + p%ModalDampingMat(k+1:k+3, :) = matmul(transpose(NodeRot), p%ModalDampingMat(k+1:k+3, :)) + p%ModalDampingMat(k+4:k+6, :) = matmul(transpose(NodeRot), p%ModalDampingMat(k+4:k+6, :)) + end do + + ! Debugging / verifying - recover the zeta values + ! phi0T_M_phi0 = matmul(transpose(eigenvectors), matmul(m%LP_MassM_LU, eigenvectors)) + ! print *, 'Verifying mass orthonormal here.' + ! phi0T_M_phi0 = matmul(transpose(eigenvectors), matmul(p%ModalDampingMat, eigenvectors)) + ! do j = 1,nDOF + ! phi0T_M_phi0(j,j) = phi0T_M_phi0(j,j) / omega(j) / 2.0_R8Ki + ! end do + + call CalcModalParticipation() + + ! Allocate memory for the velocity vector that will be multiplied by the modal damping matrix + call AllocAry(m%DampedVelocities, nDOF, 'DampedVelocities', ErrStat2, ErrMsg2); if (Failed()) return + + ! Allocate memory for the velocity vector that will be multiplied by the modal damping matrix + call AllocAry(m%ModalDampingF, nDOF, 'ModalDampingF', ErrStat2, ErrMsg2); if (Failed()) return + +contains + + logical function Failed() + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + Failed = ErrStat >= AbortErrLev + end function Failed + + subroutine CalcModalParticipation() + + REAL(R8Ki), ALLOCATABLE :: modal_participation (:, :) ! Modal participation factor + INTEGER(IntKi) :: bdModesFile ! Unit numbers for file with BD modes + + ! Theory based on Abaqus documentation + ! Only using rotational DOFs for rotations and not including + ! contributions from translations + + CALL AllocAry(modal_participation, nDOF, 6, 'modal_participation', ErrStat2, ErrMsg2) + CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'Init_ModalDamping') + + do j = 1, 6 + modal_participation(:, j) = sum(phiT_M(:, j::6)*phiT_M(:, j::6), dim=2) + end do + + ! Write to a file instead + CALL GetNewUnit(bdModesFile) + open(unit=bdModesFile, file='beamdyn_modes.csv') + + write(bdModesFile,*) '#Frequency [Hz], Zeta [Frac. Critical],'// & + 'Participation X, Participation Y, Participation Z,'// & + 'Participation RX, Participation RY, Participation RZ' + + ! Write to a file + do j = 1, numZeta + write(bdModesFile, ' (1F12.4,1F12.8,6E14.5) ') omega(j)/TwoPi_D, & + p%zeta(j), modal_participation(j, :) / sum(modal_participation(j, :)) + end do + + close(bdModesFile) + + END SUBROUTINE + +END SUBROUTINE + !----------------------------------------------------------------------------------------------------------------------------------- !> This routine is called at the end of the simulation. SUBROUTINE BD_End( u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg ) @@ -1981,7 +2191,7 @@ SUBROUTINE BD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, CALL BD_QPDataVelocity( p, x_tmp, m ) ! x%dqdt --> m%qp%vvv, m%qp%vvp ! calculate accelerations and reaction loads (in m%RHS): - CALL BD_CalcForceAcc(m%u, p, OtherState, m, ErrStat2,ErrMsg2) + CALL BD_CalcForceAcc(m%u, p, x, OtherState, m, ErrStat2,ErrMsg2) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) ELSE @@ -2105,7 +2315,7 @@ SUBROUTINE BD_CalcContStateDeriv( t, u, p, x, xd, z, OtherState, m, dxdt, ErrSta CALL BD_QPDataVelocity( p, dxdt, m ) ! x%dqdt --> m%qp%vvv, m%qp%vvp ! calculate accelerations and reaction loads (in m%RHS): - CALL BD_CalcForceAcc(m%u, p, OtherState, m, ErrStat2,ErrMsg2) + CALL BD_CalcForceAcc(m%u, p, x, OtherState, m, ErrStat2,ErrMsg2) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) if (ErrStat >= AbortErrLev) return @@ -3216,7 +3426,7 @@ SUBROUTINE BD_ElementMatrixAcc( nelem, p, OtherState, m ) CALL BD_ElasticForce( nelem, p, m, .FALSE. ) ! Calculate Fc, Fd only - IF(p%damp_flag .NE. 0) THEN + IF(p%damp_flag .EQ. 1) THEN CALL BD_DissipativeForce( nelem, p, m, .FALSE. ) ! Calculate dissipative terms on Fc, Fd ENDIF CALL BD_GravityForce( nelem, p, m, MATMUL(p%gravity,OtherState%GlbRot) ) ! Calculate Fg @@ -4955,7 +5165,7 @@ SUBROUTINE BD_ElementMatrixGA2( fact, nelem, p, OtherState, m ) CALL BD_ElasticForce( nelem,p,m,fact ) ! Calculate Fc, Fd [and if(fact): Oe, Pe, and Qe for N-R algorithm] using m%qp%E1, m%qp%RR0, m%qp%kappa, m%qp%Stif CALL BD_InertialForce( nelem,p,m,fact ) ! Calculate Fi [and Mi,Gi,Ki IF(fact)] - IF(p%damp_flag .NE. 0) THEN + IF(p%damp_flag .EQ. 1) THEN CALL BD_DissipativeForce( nelem,p,m,fact ) ! Calculate dissipative terms on Fc, Fd [and Sd, Od, Pd and Qd, betaC, Gd, Xd, Yd for N-R algorithm] ENDIF @@ -5009,7 +5219,7 @@ SUBROUTINE BD_ElementMatrixGA2( fact, nelem, p, OtherState, m ) END DO ! Dissipative terms - IF (p%damp_flag .NE. 0) THEN + IF (p%damp_flag .EQ. 1) THEN DO j=1,p%nodes_per_elem DO idx_dof2=1,p%dof_node DO i=1,p%nodes_per_elem @@ -5425,7 +5635,7 @@ SUBROUTINE BD_InitAcc( u, p, x, OtherState, m, qdotdot, ErrStat, ErrMsg ) CALL BD_QPDataVelocity(p, x, m) ! set misc vars, particularly m%RHS - CALL BD_CalcForceAcc( u, p, OtherState, m, ErrStat2, ErrMsg2 ) + CALL BD_CalcForceAcc( u, p, x, OtherState, m, ErrStat2, ErrMsg2 ) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) ! set accelerations with inputs from the root and BD_CalcForceAcc solution @@ -5509,10 +5719,11 @@ END SUBROUTINE BD_InitAcc !! !! The root reaction force is therefore calculated afterwards as !! \f$ F_\textrm{root} = f_1 - \sum_{i} m_{1,i} a_{i} \f$. -SUBROUTINE BD_CalcForceAcc( u, p, OtherState, m, ErrStat, ErrMsg ) +SUBROUTINE BD_CalcForceAcc( u, p, x, OtherState, m, ErrStat, ErrMsg ) TYPE(BD_InputType), INTENT(IN ) :: u !< Inputs at t TYPE(BD_ParameterType), INTENT(IN ) :: p !< Parameters + TYPE(BD_ContinuousStateType), INTENT(IN ) :: x !< Continuous states TYPE(BD_OtherStateType), INTENT(IN ) :: OtherState !< other states (contains ref orientation) TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables INTEGER(IntKi), INTENT( OUT) :: ErrStat !< Error status of the operation @@ -5562,6 +5773,11 @@ SUBROUTINE BD_CalcForceAcc( u, p, OtherState, m, ErrStat, ErrMsg ) ! Add force contributions from root acceleration m%LP_RHS_LU = m%LP_RHS_LU - matmul(m%LP_MassM(7:,1:6), RootAcc) + IF(p%damp_flag .EQ. 2) THEN + ! Because modal damping is already global, it wouldn't make sense in BD_AssembleRHS. + CALL BD_AddModalDampingRHS(u, p, x, OtherState, m) + ENDIF + ! Solve linear equations A * X = B for acceleration (F=ma) for nodes 2:p%node_total CALL LAPACK_getrf(n_free, n_free, m%LP_MassM_LU, m%LP_indx, ErrStat2, ErrMsg2) CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) @@ -5693,6 +5909,84 @@ SUBROUTINE BD_ComputeElementMass(nelem,p,NQPpos,EMass0_GL,elem_mass,elem_CG,elem END SUBROUTINE BD_ComputeElementMass +!----------------------------------------------------------------------------------------------------------------------------------- +!> This subroutine calculates the modal damping force +! Adds modal damping contributions to m%LP_RHS_LU +SUBROUTINE BD_AddModalDampingRHS(u, p, x, OtherState, m) + + TYPE(BD_InputType), INTENT(IN ) :: u !< Inputs at t + TYPE(BD_ParameterType), INTENT(IN ) :: p !< Parameters + TYPE(BD_ContinuousStateType), INTENT(IN ) :: x !< Continuous states + TYPE(BD_OtherStateType), INTENT(IN ) :: OtherState !< other states (contains ref orientation) + TYPE(BD_MiscVarType), INTENT(INOUT) :: m !< Misc/optimization variables + + integer(IntKi) :: j, k ! looping indexing variable + integer(IntKi) :: elem ! looping indexing for element number + integer(IntKi) :: elem_node ! looping indexing for node in the element number + real(R8Ki) :: r(3) ! nodal position relative to root + real(R8Ki) :: NodeRot(3, 3) + + ! 1. Velocities relative to root + ! element loops + do elem = 1, p%elem_total + + do elem_node = 2, p%nodes_per_elem + + ! Global node index, excluding root + j = (elem - 1) * (p%nodes_per_elem - 1) + elem_node + + ! DOF index + k = (j - 2) * 6 + + ! 1.a. Node translational velocity minus root velocity + ! x%... is at j+1 because skipping the root node and j is 1 at the first node after root. + m%DampedVelocities(k+1:k+6) = x%dqdt(:, j) - x%dqdt(:, 1) + + ! 1.b. Subtract out the rigid body rotational velocity based on the blade root rotation + + ! Vector from root to node + ! r = OtherState%GlbPos + p%uuN0(1:3, elem_node, elem) + x%q(1:3, j) - & + ! (u%RootMotion%Position(:, 1) + u%RootMotion%TranslationDisp(:, 1)) + r = p%uuN0(1:3, elem_node, elem) + x%q(1:3, j) + + m%DampedVelocities(k+1:k+3) = m%DampedVelocities(k+1:k+3) - Cross_Product(u%RootMotion%RotationVel(:, 1), r) + end do + + end do + + ! 2. rotate velocities by matmul(u%RootMotion%Orientation, OtherState%GlbRot) + ! Solve is done at the coordinate system based on time n for states at time n+1 + ! The damping matrix is technically defined for the n+1 coordinate system since that's where accel is evaluated. + ! Therefore, need to resolve the coordinate differences. + ! + ! OtherState%GlbRot = tranpose(u%RootMotion%Orientation(:, :, 1) evaluated at n) + ! here, u%RootMotion%Orientation(:, :, 1) is evaluated at n+1, but is tranposed at this point + + do j = 2, p%node_total + + ! Loop over the nodes that apply to the damping matrix, so don't include the root node. + call BD_CrvMatrixR(x%q(4:6, j), NodeRot) + + k = (j - 2) * 6 + m%DampedVelocities(k+1:k+3) = matmul(transpose(NodeRot), m%DampedVelocities(k+1:k+3)) + m%DampedVelocities(k+4:k+6) = matmul(transpose(NodeRot), m%DampedVelocities(k+4:k+6)) + end do + + ! 3. Multiply by modal damping matrix + m%ModalDampingF = matmul(p%ModalDampingMat, m%DampedVelocities) + + ! 4. Rotate to correct coordinates and subtract from m%LP_RHS_LU + do j = 2, p%node_total + + call BD_CrvMatrixR(x%q(4:6, j), NodeRot) + + k = (j - 2) * 6 + m%LP_RHS_LU(k+1:k+3) = m%LP_RHS_LU(k+1:k+3) - matmul(NodeRot, m%ModalDampingF(k+1:k+3)) + m%LP_RHS_LU(k+4:k+6) = m%LP_RHS_LU(k+4:k+6) - matmul(NodeRot, m%ModalDampingF(k+4:k+6)) + end do + +END SUBROUTINE BD_AddModalDampingRHS + subroutine BD_InitVars(u, p, x, y, m, InitOut, Linearize, ErrStat, ErrMsg) type(BD_InputType), intent(inout) :: u !< An initial guess for the input; input mesh must be defined diff --git a/modules/beamdyn/src/BeamDyn_IO.f90 b/modules/beamdyn/src/BeamDyn_IO.f90 index 5adace9939..e0322d6bb1 100644 --- a/modules/beamdyn/src/BeamDyn_IO.f90 +++ b/modules/beamdyn/src/BeamDyn_IO.f90 @@ -1056,49 +1056,28 @@ SUBROUTINE BD_ReadBladeFile(BldFile,BladeInputFileData,UnEc,ErrStat,ErrMsg) CALL ReadCom(UnIn,BldFile,'unused beam file header line 2',ErrStat2,ErrMsg2,UnEc) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - ! -------------- BLADE PARAMETER----------------------------------------------- - CALL ReadCom(UnIn,BldFile,'beam parameters',ErrStat2,ErrMsg2,UnEc) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - - CALL ReadVar(UnIn,BldFile,BladeInputFileData%station_total,'station_total','Number of blade input stations',ErrStat2,ErrMsg2,UnEc) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - - CALL AllocAry(BladeInputFileData%stiff0,6,6,BladeInputFileData%station_total,'Cross-sectional 6 by 6 stiffness matrix',ErrStat2,ErrMsg2) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL AllocAry(BladeInputFileData%mass0,6,6,BladeInputFileData%station_total,'Cross-sectional 6 by 6 mass matrix',ErrStat2,ErrMsg2) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL AllocAry(BladeInputFileData%station_eta,BladeInputFileData%station_total,'Station eta array',ErrStat2,ErrMsg2) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - - !after allocating these arrays, we'll make sure it was successful: - if (ErrStat >= AbortErrLev) then - call cleanup() - return - end if - -!FIXME: any reason not to use a logical for this? - CALL ReadVar(UnIn,BldFile,BladeInputFileData%damp_flag,'damp_flag','Damping flag',ErrStat2,ErrMsg2,UnEc) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - - ! -------------- DAMPING PARAMETER----------------------------------------------- - CALL ReadCom(UnIn,BldFile,'damping parameters',ErrStat2,ErrMsg2,UnEc) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL ReadCom(UnIn,BldFile,'mu1 to mu6',ErrStat2,ErrMsg2,UnEc) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - CALL ReadCom(UnIn,BldFile,'units',ErrStat2,ErrMsg2,UnEc) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - - CALL ReadAry(UnIn,BldFile,BladeInputFileData%beta,6,'damping coefficient','damping coefficient',ErrStat2,ErrMsg2,UnEc) - CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) - - if (ErrStat >= AbortErrLev) then - call cleanup() - return - end if - - -! -------------- DISTRIBUTED PROPERTIES-------------------------------------------- - CALL ReadCom(UnIn,BldFile,'Distributed properties',ErrStat2,ErrMsg2,UnEc) + ! Blade Parameters ---------------------------------------------------------- + CALL ReadCom(UnIn, BldFile, 'Blade Parameters', ErrStat2, ErrMsg2, UnEc); if (Failed()) return + CALL ReadVar(UnIn, BldFile, BladeInputFileData%station_total,'station_total','Number of blade input stations', ErrStat2, ErrMsg2, UnEc); if (Failed()) return + CALL AllocAry(BladeInputFileData%stiff0, 6, 6, BladeInputFileData%station_total, 'Cross-sectional 6 by 6 stiffness matrix', ErrStat2, ErrMsg2); if (Failed()) return + CALL AllocAry(BladeInputFileData%mass0, 6, 6, BladeInputFileData%station_total, 'Cross-sectional 6 by 6 mass matrix', ErrStat2, ErrMsg2); if (Failed()) return + CALL AllocAry(BladeInputFileData%station_eta, BladeInputFileData%station_total, 'Station eta array', ErrStat2, ErrMsg2); if (Failed()) return + CALL ReadVar(UnIn, BldFile, BladeInputFileData%damp_flag, 'damp_flag', 'Damping type (switch) {0: None, 1: Stiffness-Proportional, 2: Modal}', ErrStat2, ErrMsg2, UnEc); if (Failed()) return + + ! Stiffness Proportional Damping -------------------------------------------- + CALL ReadCom(UnIn, BldFile, 'Stiffness-Proportional Damping Header', ErrStat2, ErrMsg2, UnEc); if (Failed()) return + CALL ReadCom(UnIn, BldFile, 'Mu Table Header 1', ErrStat2, ErrMsg2, UnEc); if (Failed()) return + CALL ReadCom(UnIn, BldFile, 'Mu Table Header 2', ErrStat2, ErrMsg2, UnEc); if (Failed()) return + CALL ReadAry(UnIn, BldFile, BladeInputFileData%beta, 6, 'beta', 'Mu damping coefficients', ErrStat2, ErrMsg2, UnEc); if (Failed()) return + + ! Modal Damping ------------------------------------------------------------- + CALL ReadCom(UnIn, BldFile, 'Modal Damping Header', ErrStat2, ErrMsg2, UnEc); if (Failed()) return + CALL ReadVar(UnIn, BldFile, BladeInputFileData%n_modes, 'n_modes', 'Number of modal damping coefficients (-)', ErrStat2, ErrMsg2, UnEc); if (Failed()) return + CALL AllocAry(BladeInputFileData%zeta, BladeInputFileData%n_modes, 'Modal damping coefficients (zeta)', ErrStat2, ErrMsg2); if (Failed()) return + CALL ReadAry(UnIn, BldFile, BladeInputFileData%zeta, BladeInputFileData%n_modes, 'zeta', 'Damping coefficients for mode 1 through n_modes', ErrStat2, ErrMsg2, UnEc); if (Failed()) return + + ! Distributed Properties ---------------------------------------------------- + CALL ReadCom(UnIn,BldFile,'Distributed Properties',ErrStat2,ErrMsg2,UnEc) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) DO i=1,BladeInputFileData%station_total @@ -1111,7 +1090,7 @@ SUBROUTINE BD_ReadBladeFile(BldFile,BladeInputFileData,UnEc,ErrStat,ErrMsg) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) DO j=1,6 - CALL ReadAry(UnIn,BldFile,temp66(j,:),6,'siffness_matrix','Blade C/S stiffness matrix',ErrStat2,ErrMsg2,UnEc) + CALL ReadAry(UnIn,BldFile,temp66(j,:),6,'stiffness_matrix','Blade C/S stiffness matrix',ErrStat2,ErrMsg2,UnEc) CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) ENDDO if (ErrStat >= AbortErrLev) then @@ -1137,7 +1116,10 @@ SUBROUTINE BD_ReadBladeFile(BldFile,BladeInputFileData,UnEc,ErrStat,ErrMsg) return contains - !..................... + logical function Failed() + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + Failed = ErrStat >= AbortErrLev + end function subroutine cleanup() close(UnIn) return diff --git a/modules/beamdyn/src/BeamDyn_Types.f90 b/modules/beamdyn/src/BeamDyn_Types.f90 index 537c3af3a0..2122befee1 100644 --- a/modules/beamdyn/src/BeamDyn_Types.f90 +++ b/modules/beamdyn/src/BeamDyn_Types.f90 @@ -70,7 +70,9 @@ MODULE BeamDyn_Types REAL(R8Ki) , DIMENSION(:,:,:), ALLOCATABLE :: stiff0 !< C/S stiffness matrix arrays [-] REAL(R8Ki) , DIMENSION(:,:,:), ALLOCATABLE :: mass0 !< C/S mass matrix arrays [-] REAL(R8Ki) , DIMENSION(1:6) :: beta = 0.0_R8Ki !< Damping Coefficient [-] - INTEGER(IntKi) :: damp_flag = 0_IntKi !< Damping Flag: 0-No Damping, 1-Damped [-] + INTEGER(IntKi) :: n_modes = 0_IntKi !< Number of modal damping coefficients [-] + REAL(R8Ki) , DIMENSION(:), ALLOCATABLE :: zeta !< Modal damping coefficient array [-] + INTEGER(IntKi) :: damp_flag = 0_IntKi !< Damping Flag: 0-No Damping, 1-Stiffness Prop. Damped, 2-Modal Damping [-] END TYPE BladeInputData ! ======================= ! ========= BD_InputFile ======= @@ -158,6 +160,7 @@ MODULE BeamDyn_Types REAL(R8Ki) , DIMENSION(1:3) :: blade_CG = 0.0_R8Ki !< Blade center of gravity [-] REAL(R8Ki) , DIMENSION(1:3,1:3) :: blade_IN = 0.0_R8Ki !< Blade Length [-] REAL(R8Ki) , DIMENSION(1:6) :: beta = 0.0_R8Ki !< Damping Coefficient [-] + REAL(R8Ki) , DIMENSION(:), ALLOCATABLE :: zeta !< Modal Damping Coefficients [-] REAL(R8Ki) :: tol = 0.0_R8Ki !< Tolerance used in stopping criterion [-] REAL(R8Ki) , DIMENSION(:), ALLOCATABLE :: QPtN !< Quadrature (QuadPt) point locations in natural frame [-1, 1] [-] REAL(R8Ki) , DIMENSION(:), ALLOCATABLE :: QPtWeight !< Weights at each quadrature point (QuadPt) [-] @@ -212,6 +215,7 @@ MODULE BeamDyn_Types LOGICAL :: RotStates = .false. !< Orient states in rotating frame during linearization? (flag) [-] LOGICAL :: CompAeroMaps = .FALSE. !< flag to determine if BeamDyn is computing aero maps (true) or running a normal simulation (false) [-] LOGICAL :: CompAppliedLdAtRoot = .FALSE. !< flag to determine if BeamDyn should compute the applied loads at root [-] + REAL(R8Ki) , DIMENSION(:,:), ALLOCATABLE :: ModalDampingMat !< Modal damping matrix in the rotating frame [-] END TYPE BD_ParameterType ! ======================= ! ========= BD_InputType ======= @@ -305,6 +309,8 @@ MODULE BeamDyn_Types REAL(R8Ki) , DIMENSION(:), ALLOCATABLE :: LP_RHS !< Right-hand-side vector [-] REAL(R8Ki) , DIMENSION(:,:), ALLOCATABLE :: LP_StifK_LU !< Stiffness Matrix for LU [-] REAL(R8Ki) , DIMENSION(:), ALLOCATABLE :: LP_RHS_LU !< Right-hand-side vector for LU [-] + REAL(R8Ki) , DIMENSION(:), ALLOCATABLE :: DampedVelocities !< Velocity vector for applying modal damping [-] + REAL(R8Ki) , DIMENSION(:), ALLOCATABLE :: ModalDampingF !< Modal damping force in the modal damping matrix coordinates [-] INTEGER(IntKi) , DIMENSION(:), ALLOCATABLE :: LP_indx !< Index vector for LU [-] TYPE(BD_InputType) :: u !< Inputs converted to the internal BD coordinate system [-] TYPE(ModJacType) :: Jac !< Jacobian matrices and arrays corresponding to module variables [-] @@ -536,6 +542,19 @@ subroutine BD_CopyBladeInputData(SrcBladeInputDataData, DstBladeInputDataData, C DstBladeInputDataData%mass0 = SrcBladeInputDataData%mass0 end if DstBladeInputDataData%beta = SrcBladeInputDataData%beta + DstBladeInputDataData%n_modes = SrcBladeInputDataData%n_modes + if (allocated(SrcBladeInputDataData%zeta)) then + LB(1:1) = lbound(SrcBladeInputDataData%zeta) + UB(1:1) = ubound(SrcBladeInputDataData%zeta) + if (.not. allocated(DstBladeInputDataData%zeta)) then + allocate(DstBladeInputDataData%zeta(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstBladeInputDataData%zeta.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstBladeInputDataData%zeta = SrcBladeInputDataData%zeta + end if DstBladeInputDataData%damp_flag = SrcBladeInputDataData%damp_flag end subroutine @@ -555,6 +574,9 @@ subroutine BD_DestroyBladeInputData(BladeInputDataData, ErrStat, ErrMsg) if (allocated(BladeInputDataData%mass0)) then deallocate(BladeInputDataData%mass0) end if + if (allocated(BladeInputDataData%zeta)) then + deallocate(BladeInputDataData%zeta) + end if end subroutine subroutine BD_PackBladeInputData(RF, Indata) @@ -568,6 +590,8 @@ subroutine BD_PackBladeInputData(RF, Indata) call RegPackAlloc(RF, InData%stiff0) call RegPackAlloc(RF, InData%mass0) call RegPack(RF, InData%beta) + call RegPack(RF, InData%n_modes) + call RegPackAlloc(RF, InData%zeta) call RegPack(RF, InData%damp_flag) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -586,6 +610,8 @@ subroutine BD_UnPackBladeInputData(RF, OutData) call RegUnpackAlloc(RF, OutData%stiff0); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%mass0); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%beta); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%n_modes); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%zeta); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%damp_flag); if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -1197,6 +1223,18 @@ subroutine BD_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) DstParamData%blade_CG = SrcParamData%blade_CG DstParamData%blade_IN = SrcParamData%blade_IN DstParamData%beta = SrcParamData%beta + if (allocated(SrcParamData%zeta)) then + LB(1:1) = lbound(SrcParamData%zeta) + UB(1:1) = ubound(SrcParamData%zeta) + if (.not. allocated(DstParamData%zeta)) then + allocate(DstParamData%zeta(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%zeta.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%zeta = SrcParamData%zeta + end if DstParamData%tol = SrcParamData%tol if (allocated(SrcParamData%QPtN)) then LB(1:1) = lbound(SrcParamData%QPtN) @@ -1492,6 +1530,18 @@ subroutine BD_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) DstParamData%RotStates = SrcParamData%RotStates DstParamData%CompAeroMaps = SrcParamData%CompAeroMaps DstParamData%CompAppliedLdAtRoot = SrcParamData%CompAppliedLdAtRoot + if (allocated(SrcParamData%ModalDampingMat)) then + LB(1:2) = lbound(SrcParamData%ModalDampingMat) + UB(1:2) = ubound(SrcParamData%ModalDampingMat) + if (.not. allocated(DstParamData%ModalDampingMat)) then + allocate(DstParamData%ModalDampingMat(LB(1):UB(1),LB(2):UB(2)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%ModalDampingMat.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstParamData%ModalDampingMat = SrcParamData%ModalDampingMat + end if end subroutine subroutine BD_DestroyParam(ParamData, ErrStat, ErrMsg) @@ -1520,6 +1570,9 @@ subroutine BD_DestroyParam(ParamData, ErrStat, ErrMsg) if (allocated(ParamData%member_eta)) then deallocate(ParamData%member_eta) end if + if (allocated(ParamData%zeta)) then + deallocate(ParamData%zeta) + end if if (allocated(ParamData%QPtN)) then deallocate(ParamData%QPtN) end if @@ -1597,6 +1650,9 @@ subroutine BD_DestroyParam(ParamData, ErrStat, ErrMsg) if (allocated(ParamData%FEweight)) then deallocate(ParamData%FEweight) end if + if (allocated(ParamData%ModalDampingMat)) then + deallocate(ParamData%ModalDampingMat) + end if end subroutine subroutine BD_PackParam(RF, Indata) @@ -1620,6 +1676,7 @@ subroutine BD_PackParam(RF, Indata) call RegPack(RF, InData%blade_CG) call RegPack(RF, InData%blade_IN) call RegPack(RF, InData%beta) + call RegPackAlloc(RF, InData%zeta) call RegPack(RF, InData%tol) call RegPackAlloc(RF, InData%QPtN) call RegPackAlloc(RF, InData%QPtWeight) @@ -1690,6 +1747,7 @@ subroutine BD_PackParam(RF, Indata) call RegPack(RF, InData%RotStates) call RegPack(RF, InData%CompAeroMaps) call RegPack(RF, InData%CompAppliedLdAtRoot) + call RegPackAlloc(RF, InData%ModalDampingMat) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -1716,6 +1774,7 @@ subroutine BD_UnPackParam(RF, OutData) call RegUnpack(RF, OutData%blade_CG); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%blade_IN); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%beta); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%zeta); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%tol); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%QPtN); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%QPtWeight); if (RegCheckErr(RF, RoutineName)) return @@ -1794,6 +1853,7 @@ subroutine BD_UnPackParam(RF, OutData) call RegUnpack(RF, OutData%RotStates); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%CompAeroMaps); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%CompAppliedLdAtRoot); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%ModalDampingMat); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine BD_CopyInput(SrcInputData, DstInputData, CtrlCode, ErrStat, ErrMsg) @@ -2902,6 +2962,30 @@ subroutine BD_CopyMisc(SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg) end if DstMiscData%LP_RHS_LU = SrcMiscData%LP_RHS_LU end if + if (allocated(SrcMiscData%DampedVelocities)) then + LB(1:1) = lbound(SrcMiscData%DampedVelocities) + UB(1:1) = ubound(SrcMiscData%DampedVelocities) + if (.not. allocated(DstMiscData%DampedVelocities)) then + allocate(DstMiscData%DampedVelocities(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstMiscData%DampedVelocities.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstMiscData%DampedVelocities = SrcMiscData%DampedVelocities + end if + if (allocated(SrcMiscData%ModalDampingF)) then + LB(1:1) = lbound(SrcMiscData%ModalDampingF) + UB(1:1) = ubound(SrcMiscData%ModalDampingF) + if (.not. allocated(DstMiscData%ModalDampingF)) then + allocate(DstMiscData%ModalDampingF(LB(1):UB(1)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstMiscData%ModalDampingF.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstMiscData%ModalDampingF = SrcMiscData%ModalDampingF + end if if (allocated(SrcMiscData%LP_indx)) then LB(1:1) = lbound(SrcMiscData%LP_indx) UB(1:1) = ubound(SrcMiscData%LP_indx) @@ -3044,6 +3128,12 @@ subroutine BD_DestroyMisc(MiscData, ErrStat, ErrMsg) if (allocated(MiscData%LP_RHS_LU)) then deallocate(MiscData%LP_RHS_LU) end if + if (allocated(MiscData%DampedVelocities)) then + deallocate(MiscData%DampedVelocities) + end if + if (allocated(MiscData%ModalDampingF)) then + deallocate(MiscData%ModalDampingF) + end if if (allocated(MiscData%LP_indx)) then deallocate(MiscData%LP_indx) end if @@ -3103,6 +3193,8 @@ subroutine BD_PackMisc(RF, Indata) call RegPackAlloc(RF, InData%LP_RHS) call RegPackAlloc(RF, InData%LP_StifK_LU) call RegPackAlloc(RF, InData%LP_RHS_LU) + call RegPackAlloc(RF, InData%DampedVelocities) + call RegPackAlloc(RF, InData%ModalDampingF) call RegPackAlloc(RF, InData%LP_indx) call BD_PackInput(RF, InData%u) call NWTC_Library_PackModJacType(RF, InData%Jac) @@ -3158,6 +3250,8 @@ subroutine BD_UnPackMisc(RF, OutData) call RegUnpackAlloc(RF, OutData%LP_RHS); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%LP_StifK_LU); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%LP_RHS_LU); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%DampedVelocities); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%ModalDampingF); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%LP_indx); if (RegCheckErr(RF, RoutineName)) return call BD_UnpackInput(RF, OutData%u) ! u call NWTC_Library_UnpackModJacType(RF, OutData%Jac) ! Jac diff --git a/modules/beamdyn/src/Registry_BeamDyn.txt b/modules/beamdyn/src/Registry_BeamDyn.txt index 48f62ee20b..4244792344 100644 --- a/modules/beamdyn/src/Registry_BeamDyn.txt +++ b/modules/beamdyn/src/Registry_BeamDyn.txt @@ -48,12 +48,14 @@ typedef ^ BladeInputData IntKi format_index - - - #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv #the following are BDKi = R8Ki typedef ^ BladeInputData R8Ki station_eta {:} - - "Station location in eta [0,1]" -typedef ^ BladeInputData ^ stiff0 {:}{:}{:} - - "C/S stiffness matrix arrays" -typedef ^ BladeInputData ^ mass0 {:}{:}{:} - - "C/S mass matrix arrays" -typedef ^ BladeInputData ^ beta {6} - - "Damping Coefficient" - +typedef ^ BladeInputData R8Ki stiff0 {:}{:}{:} - - "C/S stiffness matrix arrays" +typedef ^ BladeInputData R8Ki mass0 {:}{:}{:} - - "C/S mass matrix arrays" +typedef ^ BladeInputData R8Ki beta {6} - - "Damping Coefficient" - +typedef ^ BladeInputData IntKi n_modes - - - "Number of modal damping coefficients" - +typedef ^ BladeInputData R8Ki zeta : - - "Modal damping coefficient array" - #end of BDKi-type variables #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -typedef ^ BladeInputData IntKi damp_flag - - - "Damping Flag: 0-No Damping, 1-Damped" +typedef ^ BladeInputData IntKi damp_flag - - - "Damping Flag: 0-No Damping, 1-Stiffness Prop. Damped, 2-Modal Damping" # ..... Input file data.............................................................................. @@ -164,6 +166,7 @@ typedef ^ ParameterType ^ blade_mass - - - typedef ^ ParameterType ^ blade_CG {3} - - "Blade center of gravity" - typedef ^ ParameterType ^ blade_IN {3}{3} - - "Blade Length" - typedef ^ ParameterType ^ beta {6} - - "Damping Coefficient" - +typedef ^ ParameterType ^ zeta {:} - - "Modal Damping Coefficients" - typedef ^ ParameterType ^ tol - - - "Tolerance used in stopping criterion" - typedef ^ ParameterType ^ QPtN {:} - - "Quadrature (QuadPt) point locations in natural frame [-1, 1]" - typedef ^ ParameterType ^ QPtWeight {:} - - "Weights at each quadrature point (QuadPt)" - @@ -222,6 +225,8 @@ typedef ^ ParameterType ^ FEweight {:}{:} typedef ^ ParameterType logical RotStates - - - "Orient states in rotating frame during linearization? (flag)" - typedef ^ ParameterType LOGICAL CompAeroMaps - .FALSE. - "flag to determine if BeamDyn is computing aero maps (true) or running a normal simulation (false)" - typedef ^ ParameterType LOGICAL CompAppliedLdAtRoot - .FALSE. - "flag to determine if BeamDyn should compute the applied loads at root" - +# .... arrays for modal damping ........................................................................................................ +typedef ^ ParameterType R8Ki ModalDampingMat {:}{:} - - "Modal damping matrix in the rotating frame" - # ..... Inputs @@ -344,6 +349,9 @@ typedef ^ MiscVarType ^ LP_MassM_LU {:}{:} - - "M typedef ^ MiscVarType ^ LP_RHS {:} - - "Right-hand-side vector" - typedef ^ MiscVarType ^ LP_StifK_LU {:}{:} - - "Stiffness Matrix for LU" - typedef ^ MiscVarType ^ LP_RHS_LU {:} - - "Right-hand-side vector for LU" - +# Velocity array for modal damping calculation +typedef ^ MiscVarType ^ DampedVelocities {:} - - "Velocity vector for applying modal damping" - +typedef ^ MiscVarType ^ ModalDampingF {:} - - "Modal damping force in the modal damping matrix coordinates" - #end of BDKi-type variables #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ typedef ^ MiscVarType IntKi LP_indx {:} - - "Index vector for LU" - diff --git a/modules/elastodyn/src/ElastoDyn.f90 b/modules/elastodyn/src/ElastoDyn.f90 index 7a0298b38e..6eb40490d2 100644 --- a/modules/elastodyn/src/ElastoDyn.f90 +++ b/modules/elastodyn/src/ElastoDyn.f90 @@ -11114,7 +11114,7 @@ subroutine ED_InitVars(u, p, x, y, m, Vars, InputFileData, Linearize, ErrStat, E do i = 1, p%NumBl call MV_AddVar(Vars%x, 'Blade'//trim(Num2LStr(i))//'Pitch', FieldAngularDisp, & DL=DatLoc(ED_x_QT), iAry=DOF_BP(i), & - Flags=VF_DerivOrder2, & + Flags=ior(VF_RotFrame, VF_DerivOrder2), & Perturb=2.0_R8Ki * D2R_D, & LinNames=['Blade pitch DOF (internal DOF index = DOF_BP('//trim(Num2LStr(i))//')), rad'], & Active=InputFileData%PitchDOF) diff --git a/modules/nwtc-library/src/NWTC_Num.f90 b/modules/nwtc-library/src/NWTC_Num.f90 index d4f36542df..a0ced39fae 100644 --- a/modules/nwtc-library/src/NWTC_Num.f90 +++ b/modules/nwtc-library/src/NWTC_Num.f90 @@ -80,7 +80,6 @@ MODULE NWTC_Num INTEGER, PARAMETER :: kernelType_TRIWEIGHT = 4 INTEGER, PARAMETER :: kernelType_TRICUBE = 5 INTEGER, PARAMETER :: kernelType_GAUSSIAN = 6 - ! constants for output formats INTEGER, PARAMETER :: Output_in_Native_Units = 0 @@ -7294,5 +7293,184 @@ SUBROUTINE fZero_R8(x, f, roots, nZeros, Period) end if END SUBROUTINE fZero_R8 +!======================================================================= + ! Copy of EigenSolve from SubDyn, migrated here to use for beamdyn modal damping. + !> Return eigenvalues, Omega, and eigenvectors + SUBROUTINE EigenSolve(K, M, N, bCheckSingularity, EigVect, Omega, ErrStat, ErrMsg ) + USE NWTC_LAPACK, only: LAPACK_ggev + + INTEGER , INTENT(IN ) :: N !< Number of degrees of freedom, size of M and K + REAL(R8Ki), INTENT(INOUT) :: K(N, N) !< Stiffness matrix + REAL(R8Ki), INTENT(INOUT) :: M(N, N) !< Mass matrix + LOGICAL, INTENT(IN ) :: bCheckSingularity ! If True, the solver will fail if rigid modes are present + REAL(R8Ki), INTENT(INOUT) :: EigVect(N, N) !< Returned Eigenvectors + REAL(R8Ki), INTENT(INOUT) :: Omega(N) !< Returned Eigenvalues + INTEGER(IntKi), INTENT( OUT) :: ErrStat !< Error status of the operation + CHARACTER(*), INTENT( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None + ! LOCALS + REAL(R8Ki), ALLOCATABLE :: WORK (:), VL(:,:), AlphaR(:), AlphaI(:), BETA(:) ! eigensolver variables + INTEGER :: i + INTEGER :: LWORK !variables for the eigensolver + INTEGER, ALLOCATABLE :: KEY(:) + INTEGER(IntKi) :: ErrStat2 + CHARACTER(ErrMsgLen) :: ErrMsg2 + REAL(R8Ki) :: normA + REAL(R8Ki) :: Omega2(N) !< Squared eigenvalues + REAL(R8Ki), parameter :: MAX_EIGENVALUE = HUGE(1.0_ReKi) ! To avoid overflow when switching to ReKi + + ErrStat = ErrID_None + ErrMsg = '' + + ! allocate working arrays and return arrays for the eigensolver + LWORK=8*N + 16 !this is what the eigensolver wants >> bjj: +16 because of MKL ?ggev documenation ( "lwork >= max(1, 8n+16) for real flavors"), though LAPACK documenation says 8n is fine + !bjj: there seems to be a memory problem in *GGEV, so I'm making the WORK array larger to see if I can figure it out + CALL AllocAry( Work, LWORK, 'Work', ErrStat2, ErrMsg2 ); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'EigenSolve') + CALL AllocAry( AlphaR, N, 'AlphaR', ErrStat2, ErrMsg2 ); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'EigenSolve') + CALL AllocAry( AlphaI, N, 'AlphaI', ErrStat2, ErrMsg2 ); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'EigenSolve') + CALL AllocAry( Beta, N, 'Beta', ErrStat2, ErrMsg2 ); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'EigenSolve') + CALL AllocAry( VL, N, N, 'VL', ErrStat2, ErrMsg2 ); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'EigenSolve') + CALL AllocAry( KEY, N, 'KEY', ErrStat2, ErrMsg2 ); if(Failed()) return + + ! --- Eigenvalue analysis + ! note: SGGEV seems to have memory issues in certain cases. The eigenvalues seem to be okay, but the eigenvectors vary wildly with different compiling options. + ! DGGEV seems to work better, so I'm making these variables R8Ki (which is set to R8Ki for now) - bjj 4/25/2014 + ! bjj: This comes from the LAPACK documentation: + ! Note: the quotients AlphaR(j)/BETA(j) and AlphaI(j)/BETA(j) may easily over- or underflow, and BETA(j) may even be zero. + ! Thus, the user should avoid naively computing the ratio Alpha/beta. However, AlphaR and AlphaI will be always less + ! than and usually comparable with norm(A) in magnitude, and BETA always less than and usually comparable with norm(B). + ! Omega2=AlphaR/BETA !Note this may not be correct if AlphaI<>0 and/or BETA=0 TO INCLUDE ERROR CHECK, also they need to be sorted + CALL LAPACK_ggev('N','V',N ,K, M, AlphaR, AlphaI, Beta, VL, EigVect, WORK, LWORK, ErrStat2, ErrMsg2) + if(Failed()) return + + ! --- Determining and sorting eigen frequencies + Omega2(:) =0.0_R8Ki + DO I=1,N !Initialize the key and calculate Omega + KEY(I)=I + !Omega2(I) = AlphaR(I)/Beta(I) + if ( EqualRealNos(real(Beta(I),ReKi),0.0_ReKi) ) then + ! --- Beta =0 + if (bCheckSingularity) call WrScr('[WARN] Large eigenvalue found, system may be ill-conditioned') + Omega2(I) = MAX_EIGENVALUE + elseif ( EqualRealNos(real(AlphaI(I),ReKi),0.0_ReKi) ) THEN + ! --- Real Eigenvalues + IF ( AlphaR(I)<0.0_R8Ki ) THEN + if ( (AlphaR(I)/Beta(I))<1e-6_R8Ki ) then + ! Tolerating very small negative eigenvalues + if (bCheckSingularity) call WrScr('[INFO] Negative eigenvalue found with small norm (system may contain rigid body mode)') + Omega2(I)=0.0_R8Ki + else + if (bCheckSingularity) call WrScr('[WARN] Negative eigenvalue found, system may be ill-conditioned.') + Omega2(I)=AlphaR(I)/Beta(I) + endif + else + Omega2(I) = AlphaR(I)/Beta(I) + endif + else + ! --- Complex Eigenvalues + normA = sqrt(AlphaR(I)**2 + AlphaI(I)**2) + if ( (normA/Beta(I))<1e-6_R8Ki ) then + ! Tolerating very small eigenvalues with imaginary part + if (bCheckSingularity) call WrScr('[WARN] Complex eigenvalue found with small norm, approximating as 0') + Omega2(I) = 0.0_R8Ki + elseif ( abs(AlphaR(I))>1e3_R8Ki*abs(AlphaI(I)) ) then + ! Tolerating very small imaginary part compared to real part... (not pretty) + if (bCheckSingularity) call WrScr('[WARN] Complex eigenvalue found with small Im compare to Re') + Omega2(I) = AlphaR(I)/Beta(I) + else + if (bCheckSingularity) call WrScr('[WARN] Complex eigenvalue found with large imaginary value)') + Omega2(I) = MAX_EIGENVALUE + endif + !call Fatal('Complex eigenvalue found, system may be ill-conditioned'); return + endif + ! Capping to avoid overflow + if (Omega2(I)> MAX_EIGENVALUE) then + Omega2(I) = MAX_EIGENVALUE + endif + enddo + + ! Sorting. LASRT has issues for double precision 64 bit on windows + !CALL ScaLAPACK_LASRT('I',N,Omega2,KEY,ErrStat2,ErrMsg2); if(Failed()) return + CALL sort_in_place(Omega2,KEY) + + ! --- Sorting eigen vectors + ! KEEP ME: scaling of the eigenvectors using generalized mass =identity criterion + ! ALLOCATE(normcoeff(N,N), STAT = ErrStat ) + ! result1 = matmul(M,EigVect) + ! result2 = matmul(transpose(EigVect),result1) + ! normcoeff=sqrt(result2) !This should be a diagonal matrix which contains the normalization factors + ! normcoeff=sqrt(matmul(transpose(EigVect),matmul(M,EigVect))) !This should be a diagonal matrix which contains the normalization factors + VL=EigVect !temporary storage for sorting EigVect + DO I=1,N + !EigVect(:,I)=VL(:,KEY(I))/normcoeff(KEY(I),KEY(I)) !reordered and normalized + EigVect(:,I)=VL(:,KEY(I)) !just reordered as Huimin had a normalization outside of this one + ENDDO + !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! + + ! --- Return Omega (capped by huge(ReKi)) and check for singularity + Omega(:) = 0.0_R8Ki + do I=1,N + if (EqualRealNos(real(Omega2(I),ReKi), 0.0_ReKi)) then ! NOTE: may be necessary for some corner numerics + Omega(i)=0.0_R8Ki + if (bCheckSingularity) then + call Fatal('Zero eigenvalue found, system may contain rigid body mode'); return + endif + elseif (Omega2(I)>0) then + Omega(i)=sqrt(Omega2(I)) + else + ! Negative eigenfrequency + print*,'>>> Wrong eigenfrequency, Omega^2=',Omega2(I) ! <<< This should never happen + Omega(i)= 0.0_R8Ki + call Fatal('Negative eigenvalue found, system may be ill-conditioned'); return + endif + enddo + + CALL CleanupEigen() + RETURN + + CONTAINS + LOGICAL FUNCTION Failed() + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'EigenSolve') + Failed = ErrStat >= AbortErrLev + if (Failed) call CleanUpEigen() + END FUNCTION Failed + + SUBROUTINE Fatal(ErrMsg_in) + character(len=*), intent(in) :: ErrMsg_in + CALL SetErrStat(ErrID_Fatal, ErrMsg_in, ErrStat, ErrMsg, 'EigenSolve'); + CALL CleanUpEigen() + END SUBROUTINE Fatal + + SUBROUTINE CleanupEigen() + IF (ALLOCATED(Work) ) DEALLOCATE(Work) + IF (ALLOCATED(AlphaR)) DEALLOCATE(AlphaR) + IF (ALLOCATED(AlphaI)) DEALLOCATE(AlphaI) + IF (ALLOCATED(Beta) ) DEALLOCATE(Beta) + IF (ALLOCATED(VL) ) DEALLOCATE(VL) + IF (ALLOCATED(KEY) ) DEALLOCATE(KEY) + END SUBROUTINE CleanupEigen + + pure subroutine sort_in_place(a,key) + real(R8Ki), intent(inout), dimension(:) :: a + integer(IntKi), intent(inout), dimension(:) :: key + integer(IntKi) :: tempI + real(R8Ki) :: temp + integer(IntKi) :: i, j + do i = 2, size(a) + j = i - 1 + temp = a(i) + tempI = key(i) + do while (j>=1 .and. a(j)>temp) + a(j+1) = a(j) + key(j+1) = key(j) + j = j - 1 + if (j<1) then + exit + endif + end do + a(j+1) = temp + key(j+1) = tempI + end do + end subroutine sort_in_place + END SUBROUTINE EigenSolve !======================================================================= END MODULE NWTC_Num diff --git a/modules/subdyn/src/FEM.f90 b/modules/subdyn/src/FEM.f90 index 84cdff9ef1..f645c10c5c 100644 --- a/modules/subdyn/src/FEM.f90 +++ b/modules/subdyn/src/FEM.f90 @@ -33,162 +33,6 @@ MODULE FEM CONTAINS -!------------------------------------------------------------------------------------------------------ -!> Return eigenvalues, Omega, and eigenvectors - -SUBROUTINE EigenSolve(K, M, N, bCheckSingularity, EigVect, Omega, ErrStat, ErrMsg ) - USE NWTC_LAPACK, only: LAPACK_ggev - INTEGER , INTENT(IN ) :: N !< Number of degrees of freedom, size of M and K - REAL(LaKi), INTENT(INOUT) :: K(N, N) !< Stiffness matrix - REAL(LaKi), INTENT(INOUT) :: M(N, N) !< Mass matrix - LOGICAL, INTENT(IN ) :: bCheckSingularity ! If True, the solver will fail if rigid modes are present - REAL(LaKi), INTENT(INOUT) :: EigVect(N, N) !< Returned Eigenvectors - REAL(LaKi), INTENT(INOUT) :: Omega(N) !< Returned Eigenvalues - INTEGER(IntKi), INTENT( OUT) :: ErrStat !< Error status of the operation - CHARACTER(*), INTENT( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None - ! LOCALS - REAL(LaKi), ALLOCATABLE :: WORK (:), VL(:,:), AlphaR(:), AlphaI(:), BETA(:) ! eigensolver variables - INTEGER :: i - INTEGER :: LWORK !variables for the eigensolver - INTEGER, ALLOCATABLE :: KEY(:) - INTEGER(IntKi) :: ErrStat2 - CHARACTER(ErrMsgLen) :: ErrMsg2 - REAL(LaKi) :: normA - REAL(LaKi) :: Omega2(N) !< Squared eigenvalues - REAL(LaKi), parameter :: MAX_EIGENVALUE = HUGE(1.0_ReKi) ! To avoid overflow when switching to ReKi - - ErrStat = ErrID_None - ErrMsg = '' - - ! allocate working arrays and return arrays for the eigensolver - LWORK=8*N + 16 !this is what the eigensolver wants >> bjj: +16 because of MKL ?ggev documenation ( "lwork >= max(1, 8n+16) for real flavors"), though LAPACK documenation says 8n is fine - !bjj: there seems to be a memory problem in *GGEV, so I'm making the WORK array larger to see if I can figure it out - CALL AllocAry( Work, LWORK, 'Work', ErrStat2, ErrMsg2 ); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'EigenSolve') - CALL AllocAry( AlphaR, N, 'AlphaR', ErrStat2, ErrMsg2 ); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'EigenSolve') - CALL AllocAry( AlphaI, N, 'AlphaI', ErrStat2, ErrMsg2 ); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'EigenSolve') - CALL AllocAry( Beta, N, 'Beta', ErrStat2, ErrMsg2 ); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'EigenSolve') - CALL AllocAry( VL, N, N, 'VL', ErrStat2, ErrMsg2 ); CALL SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,'EigenSolve') - CALL AllocAry( KEY, N, 'KEY', ErrStat2, ErrMsg2 ); if(Failed()) return - - ! --- Eigenvalue analysis - ! note: SGGEV seems to have memory issues in certain cases. The eigenvalues seem to be okay, but the eigenvectors vary wildly with different compiling options. - ! DGGEV seems to work better, so I'm making these variables LaKi (which is set to R8Ki for now) - bjj 4/25/2014 - ! bjj: This comes from the LAPACK documentation: - ! Note: the quotients AlphaR(j)/BETA(j) and AlphaI(j)/BETA(j) may easily over- or underflow, and BETA(j) may even be zero. - ! Thus, the user should avoid naively computing the ratio Alpha/beta. However, AlphaR and AlphaI will be always less - ! than and usually comparable with norm(A) in magnitude, and BETA always less than and usually comparable with norm(B). - ! Omega2=AlphaR/BETA !Note this may not be correct if AlphaI<>0 and/or BETA=0 TO INCLUDE ERROR CHECK, also they need to be sorted - CALL LAPACK_ggev('N','V',N ,K, M, AlphaR, AlphaI, Beta, VL, EigVect, WORK, LWORK, ErrStat2, ErrMsg2) - if(Failed()) return - - ! --- Determinign and sorting eigen frequencies - Omega2(:) =0.0_LaKi - DO I=1,N !Initialize the key and calculate Omega - KEY(I)=I - !Omega2(I) = AlphaR(I)/Beta(I) - if ( EqualRealNos(real(Beta(I),ReKi),0.0_ReKi) ) then - ! --- Beta =0 - if (bCheckSingularity) call WrScr('[WARN] Large eigenvalue found, system may be ill-conditioned') - Omega2(I) = MAX_EIGENVALUE - elseif ( EqualRealNos(real(AlphaI(I),ReKi),0.0_ReKi) ) THEN - ! --- Real Eigenvalues - IF ( AlphaR(I)<0.0_LaKi ) THEN - if ( (AlphaR(I)/Beta(I))<1e-6_LaKi ) then - ! Tolerating very small negative eigenvalues - if (bCheckSingularity) call WrScr('[INFO] Negative eigenvalue found with small norm (system may contain rigid body mode)') - Omega2(I)=0.0_LaKi - else - if (bCheckSingularity) call WrScr('[WARN] Negative eigenvalue found, system may be ill-conditioned.') - Omega2(I)=AlphaR(I)/Beta(I) - endif - else - Omega2(I) = AlphaR(I)/Beta(I) - endif - else - ! --- Complex Eigenvalues - normA = sqrt(AlphaR(I)**2 + AlphaI(I)**2) - if ( (normA/Beta(I))<1e-6_LaKi ) then - ! Tolerating very small eigenvalues with imaginary part - if (bCheckSingularity) call WrScr('[WARN] Complex eigenvalue found with small norm, approximating as 0') - Omega2(I) = 0.0_LaKi - elseif ( abs(AlphaR(I))>1e3_LaKi*abs(AlphaI(I)) ) then - ! Tolerating very small imaginary part compared to real part... (not pretty) - if (bCheckSingularity) call WrScr('[WARN] Complex eigenvalue found with small Im compare to Re') - Omega2(I) = AlphaR(I)/Beta(I) - else - if (bCheckSingularity) call WrScr('[WARN] Complex eigenvalue found with large imaginary value)') - Omega2(I) = MAX_EIGENVALUE - endif - !call Fatal('Complex eigenvalue found, system may be ill-conditioned'); return - endif - ! Capping to avoid overflow - if (Omega2(I)> MAX_EIGENVALUE) then - Omega2(I) = MAX_EIGENVALUE - endif - enddo - - ! Sorting. LASRT has issues for double precision 64 bit on windows - !CALL ScaLAPACK_LASRT('I',N,Omega2,KEY,ErrStat2,ErrMsg2); if(Failed()) return - CALL sort_in_place(Omega2,KEY) - - ! --- Sorting eigen vectors - ! KEEP ME: scaling of the eigenvectors using generalized mass =identity criterion - ! ALLOCATE(normcoeff(N,N), STAT = ErrStat ) - ! result1 = matmul(M,EigVect) - ! result2 = matmul(transpose(EigVect),result1) - ! normcoeff=sqrt(result2) !This should be a diagonal matrix which contains the normalization factors - ! normcoeff=sqrt(matmul(transpose(EigVect),matmul(M,EigVect))) !This should be a diagonal matrix which contains the normalization factors - VL=EigVect !temporary storage for sorting EigVect - DO I=1,N - !EigVect(:,I)=VL(:,KEY(I))/normcoeff(KEY(I),KEY(I)) !reordered and normalized - EigVect(:,I)=VL(:,KEY(I)) !just reordered as Huimin had a normalization outside of this one - ENDDO - !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++! - - ! --- Return Omega (capped by huge(ReKi)) and check for singularity - Omega(:) = 0.0_LaKi - do I=1,N - if (EqualRealNos(real(Omega2(I),ReKi), 0.0_ReKi)) then ! NOTE: may be necessary for some corner numerics - Omega(i)=0.0_LaKi - if (bCheckSingularity) then - call Fatal('Zero eigenvalue found, system may contain rigid body mode'); return - endif - elseif (Omega2(I)>0) then - Omega(i)=sqrt(Omega2(I)) - else - ! Negative eigenfrequency - print*,'>>> Wrong eigenfrequency, Omega^2=',Omega2(I) ! <<< This should never happen - Omega(i)= 0.0_LaKi - call Fatal('Negative eigenvalue found, system may be ill-conditioned'); return - endif - enddo - - CALL CleanupEigen() - RETURN - -CONTAINS - LOGICAL FUNCTION Failed() - call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, 'EigenSolve') - Failed = ErrStat >= AbortErrLev - if (Failed) call CleanUpEigen() - END FUNCTION Failed - - SUBROUTINE Fatal(ErrMsg_in) - character(len=*), intent(in) :: ErrMsg_in - CALL SetErrStat(ErrID_Fatal, ErrMsg_in, ErrStat, ErrMsg, 'EigenSolve'); - CALL CleanUpEigen() - END SUBROUTINE Fatal - - SUBROUTINE CleanupEigen() - IF (ALLOCATED(Work) ) DEALLOCATE(Work) - IF (ALLOCATED(AlphaR)) DEALLOCATE(AlphaR) - IF (ALLOCATED(AlphaI)) DEALLOCATE(AlphaI) - IF (ALLOCATED(Beta) ) DEALLOCATE(Beta) - IF (ALLOCATED(VL) ) DEALLOCATE(VL) - IF (ALLOCATED(KEY) ) DEALLOCATE(KEY) - END SUBROUTINE CleanupEigen - -END SUBROUTINE EigenSolve pure subroutine sort_in_place(a,key) real(LaKi), intent(inout), dimension(:) :: a @@ -687,6 +531,8 @@ END SUBROUTINE CraigBamptonReduction_FromPartition !! Case1: K and M are taken "as is", this is used for the "LL" part of the matrix !! Case2: K and M contain some constraints lines, and they need to be removed from the Mass/Stiffness matrix. Used for full system SUBROUTINE EigenSolveWrap(K, M, nDOF, NOmega, bCheckSingularity, EigVect, Omega, ErrStat, ErrMsg, bDOF ) + USE NWTC_Num, only: EigenSolve + INTEGER, INTENT(IN ) :: nDOF ! Total degrees of freedom of the incoming system REAL(FEKi), INTENT(IN ) :: K(nDOF, nDOF) ! stiffness matrix REAL(FEKi), INTENT(IN ) :: M(nDOF, nDOF) ! mass matrix diff --git a/openfast_io/openfast_io/FAST_reader.py b/openfast_io/openfast_io/FAST_reader.py index c7b2eb6c48..e67e7951fd 100644 --- a/openfast_io/openfast_io/FAST_reader.py +++ b/openfast_io/openfast_io/FAST_reader.py @@ -1,3 +1,4 @@ +from numbers import Number import os, re, copy import numpy as np from functools import reduce @@ -915,7 +916,7 @@ def read_BeamDynBlade(self, beamdyn_blade_file, BladeNumber = 0): f.readline() f.readline() f.readline() - #---------------------- DAMPING COEFFICIENT------------------------------------ + #---------------------- Stiffness-Proportional Damping ----------------- ln = f.readline().split() self.fst_vt['BeamDynBlade'][BladeNumber]['mu1'] = float(ln[0]) self.fst_vt['BeamDynBlade'][BladeNumber]['mu2'] = float(ln[1]) @@ -924,6 +925,11 @@ def read_BeamDynBlade(self, beamdyn_blade_file, BladeNumber = 0): self.fst_vt['BeamDynBlade'][BladeNumber]['mu5'] = float(ln[4]) self.fst_vt['BeamDynBlade'][BladeNumber]['mu6'] = float(ln[5]) f.readline() + # ------ Modal Damping [used only if damp_type=2] -------------------------------- + + self.fst_vt['BeamDynBlade'][BladeNumber]['n_modes'] = int_read(f.readline().split()[0]) + self.fst_vt['BeamDynBlade'][BladeNumber]['zeta'] = np.array([f.readline().strip().split()[:self.fst_vt['BeamDynBlade'][BladeNumber]['n_modes']]], dtype=float) + f.readline() #---------------------- DISTRIBUTED PROPERTIES--------------------------------- self.fst_vt['BeamDynBlade'][BladeNumber]['radial_stations'] = np.zeros((self.fst_vt['BeamDynBlade'][BladeNumber]['station_total'])) diff --git a/openfast_io/openfast_io/FAST_writer.py b/openfast_io/openfast_io/FAST_writer.py index 9c3eb248df..bc1bc67466 100644 --- a/openfast_io/openfast_io/FAST_writer.py +++ b/openfast_io/openfast_io/FAST_writer.py @@ -820,7 +820,7 @@ def write_BeamDynBlade(self, bldInd = None): f.write('Generated with OpenFAST_IO\n') f.write('---------------------- BLADE PARAMETERS --------------------------------------\n') f.write('{:<22} {:<11} {:}'.format(bd_blade_dict['station_total'], 'station_total', '- Number of blade input stations (-)\n')) - f.write('{:<22} {:<11} {:}'.format(bd_blade_dict['damp_type'], 'damp_type', '- Damping type: 0: no damping; 1: damped\n')) + f.write('{:<22} {:<11} {:}'.format(bd_blade_dict['damp_type'], 'damp_type', '- Damping type (switch) {0: none, 1: stiffness-proportional, 2: modal}\n')) f.write('---------------------- DAMPING COEFFICIENT------------------------------------\n') f.write(" ".join(['{:^11s}'.format(i) for i in ['mu1','mu2','mu3','mu4','mu5','mu6']])+'\n') f.write(" ".join(['{:^11s}'.format(i) for i in ['(-)','(-)','(-)','(-)','(-)','(-)']])+'\n') diff --git a/reg_tests/CTestList.cmake b/reg_tests/CTestList.cmake index 05975d2a30..c6070e21e5 100644 --- a/reg_tests/CTestList.cmake +++ b/reg_tests/CTestList.cmake @@ -358,6 +358,7 @@ of_regression("Tailfin_FreeYaw1DOF_Unsteady" "openfast;elastodyn;aerod of_regression("5MW_Land_DLL_WTurb_ADsk" "openfast;elastodyn;aerodisk") of_regression("5MW_Land_DLL_WTurb_ADsk_SED" "openfast;simple-elastodyn;aerodisk") of_regression("5MW_Land_DLL_WTurb_SED" "openfast;simple-elastodyn;aerodyn") +of_regression("IEA22MW_ModalDamping" "openfast;beamdyn;servodyn") of_aeromap_regression("5MW_Land_AeroMap" "aeromap;elastodyn;aerodyn") @@ -393,6 +394,9 @@ of_regression_linear("Fake5MW_AeroLin_B3_UA6" "-highpass=0.05" "openfas of_regression_linear("WP_Stationary_Linear" "" "openfast;linear;elastodyn") of_regression_linear("Ideal_Beam_Fixed_Free_Linear" "-highpass=0.10" "openfast;linear;beamdyn") of_regression_linear("Ideal_Beam_Free_Free_Linear" "-highpass=0.10" "openfast;linear;beamdyn") +of_regression_linear("Damped_Beam_Fixed" "-highpass=0.10" "openfast;linear;beamdyn") +of_regression_linear("Damped_Beam_Rotating" "-highpass=0.10" "openfast;linear;beamdyn") +of_regression_linear("Damped_Beam_Rotated" "-highpass=0.10" "openfast;linear;beamdyn") of_regression_linear("5MW_Land_Linear_Aero" "-highpass=0.25" "openfast;linear;elastodyn;servodyn;aerodyn") of_regression_linear("5MW_Land_Linear_Aero_CalcSteady" "-highpass=0.25" "openfast;linear;elastodyn;servodyn;aerodyn") of_regression_linear("5MW_Land_BD_Linear" "" "openfast;linear;beamdyn;servodyn") diff --git a/reg_tests/executeOpenfastLinearRegressionCase.py b/reg_tests/executeOpenfastLinearRegressionCase.py index 4af4fdf086..4f06234f92 100644 --- a/reg_tests/executeOpenfastLinearRegressionCase.py +++ b/reg_tests/executeOpenfastLinearRegressionCase.py @@ -139,7 +139,7 @@ def indent(msg, sindent='\t'): # create the local output directory if it does not already exist # and initialize it with input files for all test cases -for data in ["Ideal_Beam", "WP_Baseline"]: +for data in ["Ideal_Beam", "WP_Baseline", "Damped_Beam"]: dataDir = os.path.join(buildDirectory, data) if not os.path.isdir(dataDir): rtl.copyTree(os.path.join(moduleDirectory, data), dataDir, excludeExt=excludeExt) diff --git a/reg_tests/r-test b/reg_tests/r-test index 7d49ab184b..37cc2beb17 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 7d49ab184be35aabdbda4cbabe1d56b0fdc1e6a2 +Subproject commit 37cc2beb170b521178cff3c2e2347c0d48555329