diff --git a/docs/source/user/aerodyn/ADNodalOutputs.rst b/docs/source/user/aerodyn/ADNodalOutputs.rst index 72f7686d98..2c39eb4be0 100644 --- a/docs/source/user/aerodyn/ADNodalOutputs.rst +++ b/docs/source/user/aerodyn/ADNodalOutputs.rst @@ -29,8 +29,8 @@ channels will be named with the convention of **B**\ :math:`\mathbf{\beta}`\ the three digit node number. -Sample Nodal Outputs section -^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +**Sample Nodal Outputs section** + This sample includes the ``END`` statement from the regular outputs section. diff --git a/docs/source/user/aerodyn/examples/ad_primary_example.dat b/docs/source/user/aerodyn/examples/ad_primary_example.dat index 61e5eb13a1..8f74766f1e 100644 --- a/docs/source/user/aerodyn/examples/ad_primary_example.dat +++ b/docs/source/user/aerodyn/examples/ad_primary_example.dat @@ -1,14 +1,12 @@ -------- AERODYN v15 for OpenFAST INPUT FILE ----------------------------------------------- +------- AERODYN for OpenFAST INPUT FILE ------------------------------------------------------------- Description line that will be printed in the output file and written to the screen. ====== General Options ============================================================================ -True Echo - Echo the input to ".AD.ech"? (flag) +False Echo - Echo the input to ".AD.ech"? (flag) "default" DTAero - Time interval for aerodynamic calculations {or "default"} (s) - 1 WakeMod - Type of wake/induction model (switch) {0=none, 1=BEMT, 2=DBEMT, 3=OLAF} [WakeMod cannot be 2 or 3 when linearizing] - 1 AFAeroMod - Type of blade airfoil aerodynamics model (switch) {1=steady model, 2=Beddoes-Leishman unsteady model} [AFAeroMod must be 1 when linearizing] - 0 TwrPotent - Type of tower influence on wind based on potential flow around the tower (switch) {0=none, 1=baseline potential flow, 2=potential flow with Bak correction} + 1 Wake_Mod - Wake/induction model (switch) {0=none, 1=BEMT, 3=OLAF} [WakeMod cannot be 2 or 3 when linearizing] + 1 TwrPotent - Type tower influence on wind based on potential flow around the tower (switch) {0=none, 1=baseline potential flow, 2=potential flow with Bak correction} 0 TwrShadow - Calculate tower influence on wind based on downstream tower shadow (switch) {0=none, 1=Powles model, 2=Eames model} False TwrAero - Calculate tower aerodynamic loads? (flag) -False FrozenWake - Assume frozen wake during linearization? (flag) [used only when WakeMod=1 and when linearizing] False CavitCheck - Perform cavitation check? (flag) [AFAeroMod must be 1 when CavitCheck=true] False Buoyancy - Include buoyancy effects? (flag) False CompAA - Flag to compute AeroAcoustics calculation [used only when WakeMod = 1 or 2] @@ -19,26 +17,38 @@ False CompAA - Flag to compute AeroAcoustics calculation [us "default" SpdSound - Speed of sound in working fluid (m/s) "default" Patm - Atmospheric pressure (Pa) [used only when CavitCheck=True] "default" Pvap - Vapour pressure of working fluid (Pa) [used only when CavitCheck=True] -====== Blade-Element/Momentum Theory Options ====================================================== [unused when WakeMod=0 or 3] - 1 SkewMod - Type of skewed-wake correction model (switch) {1=uncoupled, 2=Pitt/Peters, 3=coupled} [unused when WakeMod=0 or 3] -"default" SkewModFactor - Constant used in Pitt/Peters skewed wake model {or "default" is 15/32*pi} (-) [used only when SkewMod=2; unused when WakeMod=0 or 3] -f TipLoss - Use the Prandtl tip-loss model? (flag) [unused when WakeMod=0 or 3] -f HubLoss - Use the Prandtl hub-loss model? (flag) [unused when WakeMod=0 or 3] +====== Blade-Element/Momentum Theory Options ====================================================== [unused when WakeMod=0 or 3, except for BEM_Mod] + 2 BEM_Mod - BEM model {1=legacy NoSweepPitchTwist, 2=polar} (switch) [used for all Wake_Mod to determine output coordinate system] +--- Skew correction + 1 Skew_Mod - Skew model {0=No skew model, -1=Remove non-normal component for linearization, 1=skew model active} +True SkewMomCorr - Turn the skew momentum correction on or off [used only when Skew_Mod=1] + 1 SkewRedistr_Mod - Type of skewed-wake correction model (switch) {0=no redistribution, 1=Glauert/Pitt/Peters, default=1} [used only when Skew_Mod=1] +"default" SkewRedistrFactor - Constant used in Pitt/Peters skewed wake model {or "default" is 15/32*pi} (-) [used only when Skew_Mod=1 and SkewRedistr_Mod=1] +--- BEM algorithm +True TipLoss - Use the Prandtl tip-loss model? (flag) [unused when WakeMod=0 or 3] +True HubLoss - Use the Prandtl hub-loss model? (flag) [unused when WakeMod=0 or 3] True TanInd - Include tangential induction in BEMT calculations? (flag) [unused when WakeMod=0 or 3] True AIDrag - Include the drag term in the axial-induction calculation? (flag) [unused when WakeMod=0 or 3] True TIDrag - Include the drag term in the tangential-induction calculation? (flag) [unused when WakeMod=0,3 or TanInd=FALSE] - 1E-05 IndToler - Convergence tolerance for BEMT nonlinear solve residual equation {or "default"} (-) [unused when WakeMod=0 or 3] +"Default" IndToler - Convergence tolerance for BEMT nonlinear solve residual equation {or "default"} (-) [unused when WakeMod=0 or 3] 100 MaxIter - Maximum number of iteration steps (-) [unused when WakeMod=0] -====== Dynamic Blade-Element/Momentum Theory Options ============================================== [used only when WakeMod=2] - 2 DBEMT_Mod - Type of dynamic BEMT (DBEMT) model {1=constant tau1, 2=time-dependent tau1, 3=constant tau1 with continuous formulation} (-) [used only when WakeMod=2] - 4 tau1_const - Time constant for DBEMT (s) [used only when WakeMod=2 and DBEMT_Mod=1 or 3] +--- Shear correction +False SectAvg - Use sector averaging (flag) +1 SectAvgWeighting - Weighting function for sector average {1=Uniform, default=1} within a sector centered on the blade (switch) [used only when SectAvg=True] +5 SectAvgNPoints - Number of points per sectors (-) {default=5} [used only when SectAvg=True] +-60 SectAvgPsiBwd - Backward azimuth relative to blade where the sector starts (<=0) {default=-60} (deg) [used only when SectAvg=True] + 60 SectAvgPsiFwd - Forward azimuth relative to blade where the sector ends (>=0) {default=60} (deg) [used only when SectAvg=True] +--- Dynamic wake/inflow + 2 DBEMT_Mod - Type of dynamic BEMT (DBEMT) model {0=No Dynamic Wake, -1=Frozen Wake for linearization, 1:constant tau1, 2=time-dependent tau1, 3=constant tau1 with continuous formulation} (-) + 20 tau1_const - Time constant for DBEMT (s) [used only when DBEMT_Mod=1 or 3] ====== OLAF -- cOnvecting LAgrangian Filaments (Free Vortex Wake) Theory Options ================== [used only when WakeMod=3] "unused" OLAFInputFileName - Input file for OLAF [used only when WakeMod=3] -====== Beddoes-Leishman Unsteady Airfoil Aerodynamics Options ===================================== [used only when AFAeroMod=2] - 1 UAMod - Unsteady Aero Model Switch (switch) {2=B-L Gonzalez, 3=B-L Minnema/Pierce, 4=B-L HGM 4-states, 5=B-L HGM+vortex 5 states, 6=Oye, 7=Boeing-Vertol} [used only when AFAeroMod=2] -FALSE FLookup - Flag to indicate whether a lookup for f' will be calculated (TRUE) or whether best-fit exponential equations will be used (FALSE); if FALSE S1-S4 must be provided in airfoil input files (flag) [used only when AFAeroMod=2] - 0.25 UAStartRad - Starting radius for dynamic stall (fraction of rotor radius [0.0,1.0]) [used only when AFAeroMod=2; if line is missing UAStartRad=0] - 0.95 UAEndRad - Ending radius for dynamic stall (fraction of rotor radius [0.0,1.0]) [used only when AFAeroMod=2; if line is missing UAEndRad=1] +====== Unsteady Airfoil Aerodynamics Options ==================================================== +True AoA34 - Sample the angle of attack (AoA) at the 3/4 chord or the AC point {default=True} [always used] + 3 UA_Mod - Unsteady Aero Model Switch (switch) {0=Quasi-steady (no UA), 2=B-L Gonzalez, 3=B-L Minnema/Pierce, 4=B-L HGM 4-states, 5=B-L HGM+vortex 5 states, 6=Oye, 7=Boeing-Vertol} +True FLookup - Flag to indicate whether a lookup for f' will be calculated (TRUE) or whether best-fit exponential equations will be used (FALSE); if FALSE S1-S4 must be provided in airfoil input files (flag) [used only when AFAeroMod=2] + 0 UAStartRad - Starting radius for dynamic stall (fraction of rotor radius [0.0,1.0]) [used only when AFAeroMod=2; if line is missing UAStartRad=0] + 1 UAEndRad - Ending radius for dynamic stall (fraction of rotor radius [0.0,1.0]) [used only when AFAeroMod=2; if line is missing UAEndRad=1] ====== Airfoil Information ========================================================================= 1 AFTabMod - Interpolation method for multiple airfoil tables {1=1D interpolation on AoA (first table only); 2=2D interpolation on AoA and Re; 3=2D interpolation on AoA and UserProp} (-) 1 InCol_Alfa - The column in the airfoil tables that contains the angle of attack (-) @@ -113,4 +123,4 @@ END of OutList section (the word "END" must appear in the first 3 columns of the "Vindy" "Alpha" END (the word "END" must appear in the first 3 columns of this last OutList line in the optional nodal output section) -==================================================================================================== \ No newline at end of file +==================================================================================================== diff --git a/docs/source/user/aerodyn/input.rst b/docs/source/user/aerodyn/input.rst index c5d45a51f6..55084ab942 100644 --- a/docs/source/user/aerodyn/input.rst +++ b/docs/source/user/aerodyn/input.rst @@ -3,6 +3,19 @@ Input Files =========== +Important changes introduced in v4.0 +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Some important changes have been introduced starting from version 4.0. +Please refer to :ref:`api_change_ad4x` to understand the link between the old and new inputs. + +The documentation below has been updated to incorporate these changes. + + + +Introduction +~~~~~~~~~~~~ + The user configures the aerodynamic model parameters via a primary AeroDyn input file, as well as separate input files for airfoil and blade data. When used in standalone mode, an additional driver input @@ -43,6 +56,10 @@ primary input file is given in The input file begins with two lines of header information which is for your use, but is not used by the software. + + + + General Options ~~~~~~~~~~~~~~~ @@ -65,17 +82,23 @@ for ``DTAero`` may be used to indicate that AeroDyn should employ the time step prescribed by the driver code (OpenFAST or the standalone driver program). -Set ``WakeMod`` to 0 if you want to disable rotor wake/induction effects or 1 to -include these effects using the (quasi-steady) BEM theory model. When -``WakeMod`` is set to 2, a dynamic BEM theory model (DBEMT) is used (also -referred to as dynamic inflow or dynamic wake model, see :numref:`AD_DBEMT`). When ``WakeMod`` is set -to 3, the free vortex wake model is used, also referred to as OLAF (see -:numref:`OLAF`). ``WakeMod`` cannot be set to 2 or 3 during linearization -analyses. -Set ``AFAeroMod`` to 1 to include steady blade airfoil aerodynamics or 2 -to enable UA; ``AFAeroMod`` must be 1 during linearization analyses -with AeroDyn coupled to OpenFAST. +**Wake_Mod** +Set ``Wake_Mod`` to 0 if you want to have zero induced velocities. +Set it to 1 to include these effects using the BEM theory model. +When ``Wake_Mod`` is set to 3, the free vortex wake model is used, also referred to as OLAF (see +:numref:`OLAF`). ``Wake_Mod`` cannot be set to 3 during linearization analyses. + +.. note:: + Link to old inputs: The previous input `WakeMod` is removed, `WakeMod=2` used to mean DBEMT, but this now controlled using `DBEMT_Mod`. + `WakeMod=2` is a placeholder for future induction calculation method. + + +**~~AFAeroMod~~** +This input has been removed. See ``UA_Mod`` below. + +**~~FrozenWake~~** +This input has been removed. See ``DBEMT_Mod`` below. Set ``TwrPotent`` to 0 to disable the potential-flow influence of the tower on the fluid flow local to the @@ -94,14 +117,14 @@ tower or FALSE to disable these effects. During linearization analyses with AeroDyn coupled OpenFAST and BEM enabled (``WakeMod = 1``), set the -``FrozenWake`` flag to TRUE to employ frozen-wake assumptions during -linearization (i.e. to fix the axial and tangential induces velocities, -and, at their operating-point values during linearization) or FALSE to -recalculate the induction during linearization using BEM theory. +``DBEMT_Mod=-1`` to employ frozen-wake assumptions +(i.e. to fix the axial and tangential induces velocities, and, at their operating-point values during linearization) +or +``DBEMT_Mod=3`` to use the continuous dynamic wake model. Set the ``CavitCheck`` flag to TRUE to perform a cavitation check for MHK turbines or FALSE to disable this calculation. If ``CavitCheck`` is -TRUE, ``AFAeroMod`` must be set to 1 because the cavitation check does +TRUE, ``UA_Mod`` must be set to 0 because the cavitation check does not function with unsteady airfoil aerodynamics. If ``CavitCheck`` is TRUE, the ``MHK`` flag in the AeroDyn or OpenFAST driver input file must be set to 1 or 2 to indicate an MHK turbine is being modeled. @@ -112,7 +135,7 @@ If ``Buoyancy`` is TRUE, the ``MHK`` flag in the AeroDyn or OpenFAST driver input file must be set to 1 or 2 to indicate an MHK turbine is being modeled. Set the ``CompAA`` flag to TRUE to run aero-acoustic calculations. This -option is only available for ``WakeMod = 1`` or ``2`` and is not available for +option is only available for ``Wake_Mod = 1`` and is not available for an MHK turbine. See section :numref:`AeroAcoustics` for information on how to use this feature. @@ -121,6 +144,7 @@ sub-module. See :numref:`AeroAcoustics` for information on how to use this feature. + Environmental Conditions ~~~~~~~~~~~~~~~~~~~~~~~~ @@ -147,17 +171,58 @@ around 2,000 Pa. Blade-Element/Momentum Theory Options ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -The input parameters in this section are not used when ``WakeMod = 0``. +**BEM_Mod** +Determines the kind of BEM algorithm to use. + +- ``BEM_Mod=2`` (recommended) uses the new AeroDyn BEM implementation using the local staggered polar grid coordinate system, which is more suitable for large coning. It also includes an optional momentum correction that is important for large skew (see ``SkewMomCorr``). The feature will be documented at a later time. +- ``BEM_Mod=1`` (for backward compatibility) uses the old AeroDyn BEM implementation using the NoSweepPitchTwist coordinate system. + + +.. note:: + Link to old inputs: previous implementation would have ``BEM_Mod=1`` implied. + + +.. warning:: + + ``BEM_Mod`` currently governs the coordinate system used for "ill-defined" outputs (outputs that don't have a specified coordinate system) such as the ones that ends with "x" and "y". Other ill-defined outputs are the typical BEM quantities such as "AxInd", "TnInd", "Phi", etc. These are defined in a different coordinate system depending on `BEM_Mod`. For consistency accross differents `Wake_Mod` (even when `Wake_Mod/=1`), we use `BEM_Mod` to determine the coordinate system of the ill-defined outputs. + +The following inputs in this section are only used when ``WakeMod = 1``. + + + +**Skew_Mod** +``Skew_Mod`` determines the skew correction model (for yaw and tilt): + +- ``Skew_Mod=1``: activates Glauert's skew model (recommended). This model has two components: a momentum correction (``SkewMomCorr` `), and a velocity redistribution model (``SkewRedistr_Mod``). +- ``Skew_Mod=0`` means no skew model at all (not recommended) +- ``Skew_Mod=-1`` throws away non-normal component (for linearization). This setting makes sure the wind speed is always normal to the rotor to limit periodic variation of the wind speed if the rotor is not perpendicular to the wind (e.g. tower top tilting or tilt). This is mostly needed for linearization. + +Currently (``Skew_Mod=0``) or (``Skew_Mod=1`` and ``SkewModCorr=False`` and ``SkewRedistr_Mod = 0``) are the same, both set of inputs turn off the skew correction entirely. + +.. note:: + Link to old inputs: Previous implementations always had the skew model on. `Skew_Mod=-1` replaces the old `SkewMod=0` (an option that few users were using). -``SkewMod`` determines the skewed-wake correction model. Set -``SkewMod`` to 1 to use the uncoupled BEM solution technique without -an additional skewed-wake correction. Set ``SkewMod`` to 2 to include -the Pitt/Peters correction model. **The coupled model ``SkewMod= -3`` is not available in this version of AeroDyn.** -``SkewModFactor`` is used only when ``SkewMod = 2``. Enter a scaling factor to use -in the Pitt/Peters correction model, or enter ``"default"`` to use the default -value of :math:`\frac{15 \pi}{32}`. +**SkewMomCorr** +Turns the skew momentum correction on or off [used only when ``Skew_Mod=1``] +The feature will be documented at a later time. + +.. note:: + Link to old inputs: the previous behavior would be `SkewMomCorr=False` + +**SkewRedistr_Mod** +``SkewRedistr_Mod`` allows to turn on and off the induced velocity redistribution model, and give room for other models to be selected/implemented. Default=1. + + - 0: no redistribution + - 1: Glauert (Pitt-Peters) redistribution model + +**SkewRedistrFactor** +Defines the constant used in the Glauert redistribution model (``SkewRedistr_Mod=1``). +Use ``"default"`` to use the default value of :math:`\frac{15 \pi}{32}`. + + +BEM Algorithm options +~~~~~~~~~~~~~~~~~~~~~ Set ``TipLoss`` to TRUE to include the Prandtl tip-loss model or FALSE to disable it. Likewise, set ``HubLoss`` to TRUE to include the @@ -185,22 +250,56 @@ the BEM solve is not less than or equal to ``IndToler`` in ``MaxIter``, AeroDyn will exit the BEM solver and return an error message. -Dynamic Blade-Element/Momentum Theory Options -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -The input parameters in this section are used only when ``WakeMod = 2``. +Shear corrections +~~~~~~~~~~~~~~~~~ + +The BEM algorithm may need to be corrected to account for shear. +Currently, a sector average correction is implemented, as a beta feature, to limit fluctuations associated with variations of wind speed as the blade rotates. + +The feature will be documented at a later time and is still at an experimental stage. + +**SectAvg** Use Sector Averaging (flag). +The method uses sectors expanding forward and backward relative to the current azimuth of the blade (see ``SectAvgPsiBwd`` and ``SectAvgPsiFwd``). +The velocity is averaged within this sector by attributing different weighting at different points in the sector (see ``SectAvgWeighting``). + +**SectAvgWeighting** Weighting function for sector average. +1=Uniform (switch) [used only when ``SectAvg=True``]. Default is 1. + +**SectAvgNPoints** Number of points per sectors (-) [used only when ``SectAvg=True``]. Default is 5. + +**SectAvgPsiBwd** Backward azimuth (in degrees) relative to the blade azimuth where the sector starts. Must be negative. [used only when SectAvg=True]. Default is -60 deg. + +**SectAvgPsiFwd** Forward azimuth (in degrees) relative to the blade azimuth where the sector ends. Must be positive. [used only when SectAvg=True]. Default is 60 deg. + + + + + +Dynamic Wake / Dynamic inflow model +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The input parameters in this section are used only when ``Wake_Mod = 1``. The theory is described in :numref:`AD_DBEMT`. -There are three options available for ``DBEMT_Mod``: +The dynamic wake (also called dynamic inflow) model is governed by the input ``DBEMT_Mod``: +- ``0``: no dynamic wake, also called quasi-steady wake model (not recommended). +- ``-1``: frozen wake, the induced velocities at a given operating point will remain constant (useful for simplified linearization only). - ``1``: discrete-time Oye's model, with constant :math:`\tau_1` - ``2``: discrete-time Oye's model, with varying :math:`\tau_1`, automatically adjusted based on inflow. (recommended for time-domain simulations) - ``3``: continuous-time Oye's model, with constant :math:`\tau_1` (recommended for linearization) For ``DBEMT_Mod=1`` or ``DBEMT_Mod=3`` it is the user responsability to set the value of :math:`\tau_1` (i.e. ``tau1_const``) according to the expression given in :numref:`AD_DBEMT`, using an estimate of what the mean axial induction (:math:`\overline{a}`) and the mean relative wind velocity across the rotor (:math:`\overline{U_0}`) are for a given simulation. -The option ``DBEMT_Mod=3`` is the only one that can be used for linearization. +Only the options ``DBEMT_Mod={-1,3}`` can be used for linearization. + +.. note:: + Link to old inputs: + The option `DBEMT_Mod=-1` has the same behavior as the old `FrozenWake=True`. + `DBEMT_Mod=0` has the same behavior as the previous `WakeMod=1` option. + `DBEMT_Mod=J` (`J` in `1,2,3`) , has the same behavior as the previous `WakeMod=2` & `DBEMT_Mod=J` @@ -219,11 +318,23 @@ for this input file. Unsteady Airfoil Aerodynamics Options ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -The input parameters in this section are used only when ``AFAeroMod -= 2``. -``UAMod`` determines the UA model. It has the following options: +**AoA34** +Determine whether the baseline angle of attack is sampled at the 3/4 chord or at the aerodynamic center point. +Most ``UA_Mod`` will require `AoA34` to be set to true. But when using quasi-steady aerodynamics, the user may want to set it to true or false. + +.. warning:: + This feature is currently not implemented due to a lag between the `dev` and `dev-unstable` branch. + +.. note:: + Link to previous inputs: `AFAeroMod=1` implies `AoA34=False`. But to have a fair comparison between steady and unsteady aerodynamics model, it would be best to set `AoA34=True` when doing quasi-steady aero. + + + +``UA_Mod`` determines the UA model. It has the following options: + +- ``0``: no unsteady arifoil aerodynamics, - ``1``: the discrete-time model of Beddoes-Leishman (B-L) (**not currently functional**), - ``2``: the extensions to B-L developed by González (changes in Cn, Cc, Cm) - ``3``: the extensions to B-L developed by Minnema/Pierce (changes in Cc and Cm) @@ -232,8 +343,10 @@ The input parameters in this section are used only when ``AFAeroMod - ``6``: 1-state continuous-time developed by Oye - ``7``: discrete-time Boeing-Vertol (BV) model -Linearization is supported with ``UAMod=4,5,6`` (which use continuous-time states) but not with the other models. The different models are described in :numref:`AD_UA`. +Linearization is supported with ``UA_Mod=4,5,6`` (which use continuous-time states) but not with the other models. The different models are described in :numref:`AD_UA`. +.. note:: + Link to old inputs: If `UA_Mod>0`, then this is equivalent to the old `AFAeroMod=2`. **While all of the UA models are documented in this manual, the original B-L model is not yet functional. Testing has shown @@ -256,7 +369,7 @@ is determined via a lookup into the static lift-force coefficient and drag-force coefficient data. **Using best-fit exponential equations (``FLookup = FALSE``) is not yet available, so ``FLookup`` must be ``TRUE`` in this version of AeroDyn.** Note, ``FLookup`` is not used -when ``UAMod=4`` or ``UAMod=5``. +when ``UA_Mod=4`` or ``UA_Mod=5``. ``UAStartRad`` is the starting rotor radius where dynamic stall will be turned on. Enter a number between 0 and 1, representing a fraction of rotor radius, @@ -528,7 +641,7 @@ if the keyword ``DEFAULT`` is entered in place of a numerical value, ``RelThickness`` is the non-dimensional thickness of the airfoil (thickness over chord ratio), expressed as a fraction (not a percentage), typically between 0.1 and 1. -The parameter is currently used when ``UAMod=7``, but might be used more in the future. +The parameter is currently used when ``UA_Mod=7``, but might be used more in the future. The default value of 0.2 if provided for convenience. ``NonDimArea`` is the nondimensional airfoil area (normalized by the @@ -583,24 +696,24 @@ or calculating it based on the polar coefficient data in the airfoil table: - ``alphaUpper`` specifies the AoA (in degrees) of the upper boundary of fully-attached region of the cn or cl curve. It is used to - compute the separation function when ``UAMod=5``. + compute the separation function when ``UA_Mod=5``. - ``alphaLower`` specifies the AoA (in degrees) of the lower boundary of fully-attached region of the cn or cl curve. It is used to - compute the separation function when ``UAMod=5``. (The separation function + compute the separation function when ``UA_Mod=5``. (The separation function will have a value of 1 between ``alphaUpper`` and ``alphaLower``.) - ``eta_e`` is the recovery factor and typically has a value in the - range [0.85 to 0.95] for ``UAMod = 1``; if the keyword ``DEFAULT`` is + range [0.85 to 0.95] for ``UA_Mod = 1``; if the keyword ``DEFAULT`` is entered in place of a numerical value, ``eta_e`` is set to 0.9 for - ``UAMod = 1``, but ``eta_e`` is set to 1.0 for other ``UAMod`` + ``UA_Mod = 1``, but ``eta_e`` is set to 1.0 for other ``UA_Mod`` values and whenever ``FLookup = TRUE``; - ``C_nalpha`` is the slope of the 2D normal force coefficient curve in the linear region; - ``C_lalpha`` is the slope of the 2D normal lift coefficient curve - in the linear region; Used for ``UAMod=4,6``. + in the linear region; Used for ``UA_Mod=4,6``. - ``T_f0`` is the initial value of the time constant associated with *Df* in the expressions of *Df* and *f’*; if the keyword ``DEFAULT`` is @@ -661,19 +774,19 @@ or calculating it based on the polar coefficient data in the airfoil table: to 1, based on experimental results; - ``S1`` is the constant in the best fit curve of *f* for - ``alpha0`` :math:`\le` AoA :math:`\le` ``alpha1`` for ``UAMod = 1`` (and is unused + ``alpha0`` :math:`\le` AoA :math:`\le` ``alpha1`` for ``UA_Mod = 1`` (and is unused otherwise); by definition, it depends on the airfoil; - ``S2`` is the constant in the best fit curve of *f* for AoA > - ``alpha1`` for ``UAMod = 1`` (and is unused otherwise); by + ``alpha1`` for ``UA_Mod = 1`` (and is unused otherwise); by definition, it depends on the airfoil; - ``S3`` is the constant in the best fit curve of *f* for - ``alpha2`` :math:`\le` AoA :math:`\le` ``alpha0`` for ``UAMod = 1`` (and is unused + ``alpha2`` :math:`\le` AoA :math:`\le` ``alpha0`` for ``UA_Mod = 1`` (and is unused otherwise); by definition, it depends on the airfoil; - ``S4`` is the constant in the best fit curve of *f* for AoA < - ``alpha2`` for ``UAMod = 1`` (and is unused otherwise); by + ``alpha2`` for ``UA_Mod = 1`` (and is unused otherwise); by definition, it depends on the airfoil; - ``Cn1`` is the critical value of :math:`C^{\prime}_n` at leading-edge separation for @@ -700,22 +813,22 @@ or calculating it based on the polar coefficient data in the airfoil table: location at zero-lift AoA, positive for nose up; - ``k0`` is a constant in the best fit curve of :math:`\hat{x}_{cp}` and equals for :math:`\hat{x}_{AC}-0.25` - ``UAMod = 1`` (and is unused otherwise); + ``UA_Mod = 1`` (and is unused otherwise); -- ``k1`` is a constant in the best fit curve of :math:`\hat{x}_{cp}` for ``UAMod = 1`` +- ``k1`` is a constant in the best fit curve of :math:`\hat{x}_{cp}` for ``UA_Mod = 1`` (and is unused otherwise); -- ``k2`` is a constant in the best fit curve of :math:`\hat{x}_{cp}` for ``UAMod = 1`` +- ``k2`` is a constant in the best fit curve of :math:`\hat{x}_{cp}` for ``UA_Mod = 1`` (and is unused otherwise); -- ``k3`` is a constant in the best fit curve of :math:`\hat{x}_{cp}` for ``UAMod = 1`` +- ``k3`` is a constant in the best fit curve of :math:`\hat{x}_{cp}` for ``UA_Mod = 1`` (and is unused otherwise); - ``k1_hat`` is a constant in the expression of *Cc* due to - leading-edge vortex effects for ``UAMod = 1`` (and is unused + leading-edge vortex effects for ``UA_Mod = 1`` (and is unused otherwise); -- ``x_cp_bar`` is a constant in the expression of :math:`\hat{x}_{cp}^{\nu}` for ``UAMod = +- ``x_cp_bar`` is a constant in the expression of :math:`\hat{x}_{cp}^{\nu}` for ``UA_Mod = 1`` (and is unused otherwise); if the keyword ``DEFAULT`` is entered in place of a numerical value, ``x_cp_bar`` is set to 0.2; and diff --git a/docs/source/user/aerodyn/modeling.rst b/docs/source/user/aerodyn/modeling.rst index fa245af95c..6451958cfc 100644 --- a/docs/source/user/aerodyn/modeling.rst +++ b/docs/source/user/aerodyn/modeling.rst @@ -110,10 +110,12 @@ Linearization When coupled to FAST, AeroDyn can be linearized as part of the -linearization of the full coupled solution. When induction is enabled -(``WakeMod = 1``), we recommend to base the linearized solution on the -frozen-wake assumption, by setting ``FrozenWake = TRUE``. The UA -models are not set up to support linearization, so, UA must be disabled -during linearization by setting ``AFAeroMod = 1``. Linearization is not -currently possible when modeling an MHK turbine, but we will attempt to -enable it in an upcoming release. +linearization of the full coupled solution. +A subset of the AeroDyn modules options are available. + +Dynamic wake can be linearized with +`DBEMT_Mod=-1` (frozen-wake) +`DBEMT_Mod=3` (dynamic continuous state-space model). + +Unsteady aerodynamics can be linearized with: +`UAMod={0, 4, 5, 7`} (no UA, or continuous state-space models). diff --git a/docs/source/user/api_change.rst b/docs/source/user/api_change.rst index 5cf14c37a2..bc34500391 100644 --- a/docs/source/user/api_change.rst +++ b/docs/source/user/api_change.rst @@ -33,6 +33,44 @@ SubDyn 59\* \*Exact line number depends on number of entries in various preceeding tables. +.. _api_change_ad4x: + +AeroDyn changes starting from v4.x +---------------------------------- + +The table below shows how to convert from the Old AeroDyn inputs to the new AeroDyn inputs. +Additional ressources: + +- The AeroDyn input file description (:numref:`ad_input`) for more details on the new inputs. + +- The `discussion `__ that led to these new inputs. + +- An example of AeroDyn input file at it's latest format: :download:`Example `: + +- A directory with a working example: `here `__ + + +=========================== ========================================================= +Old inputs Corresponding new inputs +=========================== ========================================================= +`WakeMod=0` `Wake_Mod=0` +`WakeMod=1` ("BEM") `Wake_Mod=1` and `DBEMT_Mod=0` and `BEM_Mod=1` +`WakeMod=2` ("DBEMT") `Wake_Mod=1` and `DBEMT_Mod={1,2,3}` +`WakeMod=3` ("OLAF") `Wake_Mod=3` +`AFAeroMod=1` `UA_Mod=0` and `AoA34=False` +`AFAeroMod=2` `UA_Mod>0` and `AoA34=False` and `UA_Mod=UAMod` +`FrozenWake=True` `DBEMT_Mod=-1` +`FrozenWake=False` `DBEMT_Mod=0` (quasi-steady) or `DBEMT_Mod>0` (dynamic) +`SkewMod=2` (Glauert) `Skew_Mod=1` and `SkewRedistr_Mod=1` +`SkewMod=0` (Orthogonal) `Skew_Mod=-1` +`SkewModFactor` `SkewRedistrFactor` +`UAMod={2-7}` `UA_Mod={2-7}` and `AoA34=True` +=========================== ========================================================= + + + + + OpenFAST v3.5.2 to OpenFAST v3.5.3 ---------------------------------- diff --git a/modules/aerodyn/src/AeroDyn.f90 b/modules/aerodyn/src/AeroDyn.f90 index 81669a3056..83fd79f257 100644 --- a/modules/aerodyn/src/AeroDyn.f90 +++ b/modules/aerodyn/src/AeroDyn.f90 @@ -318,20 +318,22 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut CALL ParsePrimaryFileInfo( PriPath, InitInp, InitInp%InputFile, p%RootName, NumBlades, interval, FileInfo_In, InputFileData, UnEcho, ErrStat2, ErrMsg2 ) if (Failed()) return; - ! Temporary HACK, for WakeMod=10, 11 or 12 use AeroProjMod 2 (will trigger PolarBEM) - if (InputFileData%WakeMod==10) then - call WrScr(' WARNING: WakeMod=10 is a temporary hack. Using new projection method with WakeMod=0.') - InputFileData%WakeMod = 0 - AeroProjMod(:) = 2 - elseif (InputFileData%WakeMod==11) then - call WrScr(' WARNING: WakeMod=11 is a temporary hack. Using new projection method with WakeMod=1.') - InputFileData%WakeMod = 1 - AeroProjMod(:) = 2 - elseif (InputFileData%WakeMod==12) then - call WrScr(' WARNING: WakeMod=12 is a temporary hack. Using new projection method with WakeMod=2.') - InputFileData%WakeMod = 2 - AeroProjMod(:) = 2 - endif + ! --- "Automatic handling of AeroProjMod + do iR = 1, nRotors + if (AeroProjMod(iR) == -1) then + if (InputFileData%Wake_Mod /= WakeMod_BEMT) then + ! For BEMT, we don't throw a warning + call WrScr('[INFO] Using the input file input `BEM_Mod` to match BEM coordinate system outputs') + endif + select case (InputFileData%BEM_Mod) + case (BEMMod_2D); AeroProjMod(ir) = APM_BEM_NoSweepPitchTwist + case (BEMMod_3D); AeroProjMod(ir) = APM_BEM_Polar + case default; call Fatal('Input `BEM_Mod` not supported: '//trim(num2lstr(InputFileData%BEM_Mod))); return + end select + + endif + enddo + ! ----------------------------------------------------------------- ! Read the AeroDyn blade files, or copy from passed input @@ -343,20 +345,13 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut ! bjj: do we put a warning here if any of these values aren't currently set this way? if (InitInp%CompAeroMaps) then InputFileData%DTAero = interval ! we're not using this, so set it to something "safe" - do iR = 1, nRotors - InputFileData%AFAeroMod = AFAeroMod_Steady - InputFileData%TwrPotent = TwrPotent_none - InputFileData%TwrShadow = TwrShadow_none - InputFileData%TwrAero = .false. - InputFileData%FrozenWake = .false. - !InputFileData%CavitCheck = .false. - !InputFileData%TFinAero = .false. ! not sure if this needs to be set or not - end do - - if (InputFileData%WakeMod == WakeMod_DBEMT) then - ! these models (DBEMT and BEMT) should be the same at the first time step, so we'll simplify here - InputFileData%WakeMod = WakeMod_BEMT - end if + InputFileData%UA_Mod = UA_None + InputFileData%TwrPotent = TwrPotent_none + InputFileData%TwrShadow = TwrShadow_none + InputFileData%TwrAero = .false. + !InputFileData%CavitCheck = .false. + !InputFileData%TFinAero = .false. ! not sure if this needs to be set or not + InputFileData%DBEMT_Mod = DBEMT_none end if ! Validate the inputs @@ -373,11 +368,10 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut ! set the rest of the parameters - p%SkewMod = InputFileData%SkewMod + p%Skew_Mod = InputFileData%Skew_Mod do iR = 1, nRotors - !p%rotors(iR)%AeroProjMod = InitInp%rotors(iR)%AeroProjMod p%rotors(iR)%AeroProjMod = AeroProjMod(iR) - p%rotors(iR)%AeroBEM_Mod = InitInp%rotors(iR)%AeroBEM_Mod + call WrScr(' AeroDyn: projMod: '//trim(num2lstr(p%rotors(iR)%AeroProjMod))) call SetParameters( InitInp, InputFileData, InputFileData%rotors(iR), p%rotors(iR), p, ErrStat2, ErrMsg2 ) if (Failed()) return; enddo @@ -422,7 +416,7 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut ! initialize BEMT after setting parameters and inputs because we are going to use the already- ! calculated node positions from the input meshes - if (p%WakeMod /= WakeMod_FVW) then + if (p%Wake_Mod /= WakeMod_FVW) then do iR = 1, nRotors call Init_BEMTmodule( InputFileData, InputFileData%rotors(iR), u%rotors(iR), m%rotors(iR)%BEMT_u(1), p%rotors(iR), p, x%rotors(iR)%BEMT, xd%rotors(iR)%BEMT, z%rotors(iR)%BEMT, & OtherState%rotors(iR)%BEMT, m%rotors(iR)%BEMT_y, m%rotors(iR)%BEMT, ErrStat2, ErrMsg2 ) @@ -441,7 +435,7 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut end if enddo - else ! if (p%WakeMod == WakeMod_FVW) then + else ! if (p%Wake_Mod == WakeMod_FVW) then !------------------------------------------------------------------------------------------------- ! Initialize FVW module if it is used @@ -477,7 +471,7 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut ! many states are in the BEMT module, which were initialized in BEMT_Init() do iR = 1, nRotors - call Init_MiscVars(m%rotors(iR), p%rotors(iR), u%rotors(iR), y%rotors(iR), errStat2, errMsg2) + call Init_MiscVars(m%rotors(iR), p%rotors(iR), p, u%rotors(iR), y%rotors(iR), errStat2, errMsg2) if (Failed()) return; enddo @@ -549,6 +543,12 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut call Cleanup() contains + subroutine Fatal(errMsg_in) + character(*), intent(in) :: errMsg_in + call SetErrStat(ErrID_Fatal, errMsg_in, ErrStat, ErrMsg, RoutineName ) + call Cleanup() + end subroutine Fatal + logical function Failed() CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName ) Failed = ErrStat >= AbortErrLev @@ -599,7 +599,7 @@ subroutine AD_ReInit(p, x, xd, z, OtherState, m, Interval, ErrStat, ErrMsg ) ! and the UA filter end if - if (p%WakeMod /= WakeMod_FVW) then + if (p%Wake_Mod /= WakeMod_FVW) then do IR=1, size(p%rotors) call BEMT_ReInit(p%rotors(iR)%BEMT,x%rotors(iR)%BEMT,xd%rotors(iR)%BEMT,z%rotors(iR)%BEMT,OtherState%rotors(iR)%BEMT,m%rotors(iR)%BEMT,ErrStat,ErrMsg) @@ -617,9 +617,10 @@ subroutine AD_ReInit(p, x, xd, z, OtherState, m, Interval, ErrStat, ErrMsg ) end subroutine AD_ReInit !---------------------------------------------------------------------------------------------------------------------------------- !> This routine initializes (allocates) the misc variables for use during the simulation. -subroutine Init_MiscVars(m, p, u, y, errStat, errMsg) +subroutine Init_MiscVars(m, p, p_AD, u, y, errStat, errMsg) type(RotMiscVarType), intent(inout) :: m !< misc/optimization data (not defined in submodules) type(RotParameterType), intent(in ) :: p !< Parameters + type(AD_ParameterType), intent(in ) :: p_AD !< Parameters type(RotInputType), intent(inout) :: u !< input for HubMotion mesh (create sibling mesh here) type(RotOutputType), intent(inout) :: y !< output (create mapping between output and otherstate mesh here) integer(IntKi), intent( out) :: errStat !< Error status of the operation @@ -639,6 +640,9 @@ subroutine Init_MiscVars(m, p, u, y, errStat, errMsg) call AllocAry( m%DisturbedInflow, 3_IntKi, p%NumBlNds, p%numBlades, 'm%DisturbedInflow', ErrStat2, ErrMsg2 ) ! must be same size as u%InflowOnBlade call SetErrStat( errStat2, errMsg2, errStat, errMsg, RoutineName ) +if ((p_AD%SectAvg) .and. ((p_AD%Wake_Mod == WakeMod_BEMT)) ) then + call AllocAry( m%SectAvgInflow, 3_IntKi, p%NumBlNds, p%numBlades, 'm%SectAvgInflow' , ErrStat2, ErrMsg2 ); if(Failed()) return +endif call AllocAry( m%orientationAnnulus, 3_IntKi, 3_IntKi, p%NumBlNds, p%numBlades, 'm%orientationAnnulus', ErrStat2, ErrMsg2 ) call SetErrStat( errStat2, errMsg2, errStat, errMsg, RoutineName ) call AllocAry( m%R_li, 3_IntKi, 3_IntKi, p%NumBlNds, p%numBlades, 'm%R_li', ErrStat2, ErrMsg2 ) @@ -903,6 +907,12 @@ subroutine Init_MiscVars(m, p, u, y, errStat, errMsg) end if m%FirstWarn_TowerStrike = .true. + +contains + logical function Failed() + call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + Failed = ErrStat >= AbortErrLev + end function Failed end subroutine Init_MiscVars !---------------------------------------------------------------------------------------------------------------------------------- @@ -1326,25 +1336,26 @@ subroutine SetParameters( InitInp, InputFileData, RotData, p, p_AD, ErrStat, Err ! NOTE: p_AD%FlowField is set in the glue code (or ADI module); seems like FlowField should be an initialization input so that would be clearer for new developers... - p_AD%UA_Flag = InputFileData%AFAeroMod == AFAeroMod_BL_unsteady + p_AD%UA_Flag = InputFileData%UA_Mod > UA_None p_AD%CompAeroMaps = InitInp%CompAeroMaps + p_AD%SectAvg = InputFileData%SectAvg + p_AD%SA_Weighting = InputFileData%SA_Weighting + p_AD%SA_PsiBwd = InputFileData%SA_PsiBwd*D2R + p_AD%SA_PsiFwd = InputFileData%SA_PsiFwd*D2R + p_AD%SA_nPerSec = InputFileData%SA_nPerSec + p%MHK = InitInp%MHK p_AD%DT = InputFileData%DTAero - p_AD%WakeMod = InputFileData%WakeMod + p_AD%Wake_Mod = InputFileData%Wake_Mod + p%DBEMT_Mod = InputFileData%DBEMT_Mod p%TwrPotent = InputFileData%TwrPotent p%TwrShadow = InputFileData%TwrShadow p%TwrAero = InputFileData%TwrAero p%CavitCheck = InputFileData%CavitCheck p%Buoyancy = InputFileData%Buoyancy - - if (InitInp%Linearize .and. InputFileData%WakeMod == WakeMod_BEMT) then - p%FrozenWake = InputFileData%FrozenWake - else - p%FrozenWake = .FALSE. - end if p%CompAA = InputFileData%CompAA @@ -1563,7 +1574,7 @@ subroutine AD_End( u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg ) ! Place any last minute operations or calculations here: ! End the FVW submodule - if (p%WakeMod == WakeMod_FVW ) then + if (p%Wake_Mod == WakeMod_FVW ) then if ( p%UA_Flag ) then do iW=1,p%FVW%nWings @@ -1661,14 +1672,14 @@ subroutine AD_UpdateStates( t, n, u, utimes, p, x, xd, z, OtherState, m, errStat call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) do iR = 1,size(p%rotors) - call SetInputs(p%rotors(iR), p, uInterp%rotors(iR), m%rotors(iR), i, errStat2, errMsg2) + call SetInputs(t, p%rotors(iR), p, uInterp%rotors(iR), m%rotors(iR), i, errStat2, errMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) enddo enddo - if (p%WakeMod /= WakeMod_FVW) then + if (p%Wake_Mod /= WakeMod_FVW) then do iR = 1,size(p%rotors) ! Call into the BEMT update states NOTE: This is a non-standard framework interface!!!!! GJH call BEMT_UpdateStates(t, n, m%rotors(iR)%BEMT_u(:), BEMT_utimes, p%rotors(iR)%BEMT, x%rotors(iR)%BEMT, xd%rotors(iR)%BEMT, z%rotors(iR)%BEMT, OtherState%rotors(iR)%BEMT, p%AFI, m%rotors(iR)%BEMT, errStat2, errMsg2) @@ -1872,7 +1883,7 @@ subroutine AD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg, if(Failed()) return enddo - if (p%WakeMod == WakeMod_FVW) then + if (p%Wake_Mod == WakeMod_FVW) then ! This needs to extract the inputs from the AD data types (mesh) and copy pieces for the FVW module call SetInputsForFVW(p, (/u/), m, errStat2, errMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) @@ -1954,10 +1965,10 @@ subroutine RotCalcOutput( t, u, p, p_AD, x, xd, z, OtherState, y, m, m_AD, iRot, CalcWriteOutput = .true. ! by default, calculate WriteOutput unless told that we do not need it end if - call SetInputs(p, p_AD, u, m, indx, errStat2, errMsg2) + call SetInputs(t, p, p_AD, u, m, indx, errStat2, errMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then ! Call the BEMT module CalcOutput. Notice that the BEMT outputs are purposely attached to AeroDyn's MiscVar structure to ! avoid issues with the coupling code @@ -2089,10 +2100,10 @@ subroutine AD_CavtCrit(u, p, m, errStat, errMsg) do j = 1,p%rotors(iR)%numBlades ! Loop through all blades do i = 1,p%rotors(iR)%NumBlNds ! Loop through all nodes - if ( p%WakeMod == WakeMod_BEMT .or. p%WakeMod == WakeMod_DBEMT ) then + if ( p%Wake_Mod == WakeMod_BEMT ) then Vreltemp = m%rotors(iR)%BEMT_y%Vrel(i,j) Cpmintemp = m%rotors(iR)%BEMT_y%Cpmin(i,j) - else if ( p%WakeMod == WakeMod_FVW ) then + else if ( p%Wake_Mod == WakeMod_FVW ) then iW = p%FVW%Bld2Wings(iR,j) Vreltemp = m%FVW%W(iW)%BN_Vrel(i) Cpmintemp = m%FVW%W(iW)%BN_Cpmin(i) @@ -2581,7 +2592,7 @@ subroutine RotCalcConstrStateResidual( Time, u, p, p_AD, x, xd, z, OtherState, m end if - call SetInputs(p, p_AD, u, m, indx, errStat2, errMsg2) + call SetInputs(Time, p, p_AD, u, m, indx, errStat2, errMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) @@ -2621,7 +2632,7 @@ subroutine RotCalcContStateDeriv( t, u, p, p_AD, x, xd, z, OtherState, m, dxdt, ErrStat = ErrID_None ErrMsg = "" - call SetInputs(p, p_AD, u, m, InputIndex, ErrStat2, ErrMsg2) + call SetInputs(t, p, p_AD, u, m, InputIndex, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) call BEMT_CalcContStateDeriv( t, m%BEMT_u(InputIndex), p%BEMT, x%BEMT, xd%BEMT, z%BEMT, OtherState%BEMT, m%BEMT, dxdt%BEMT, p_AD%AFI, ErrStat2, ErrMsg2 ) @@ -2631,7 +2642,8 @@ END SUBROUTINE RotCalcContStateDeriv !---------------------------------------------------------------------------------------------------------------------------------- !> This subroutine converts the AeroDyn inputs into values that can be used for its submodules. It calculates the disturbed inflow !! on the blade if tower shadow or tower influence are enabled, then uses these values to set m%BEMT_u(indx). -subroutine SetInputs(p, p_AD, u, m, indx, errStat, errMsg) +subroutine SetInputs(t, p, p_AD, u, m, indx, errStat, errMsg) + real(DbKi), intent(in ) :: t !< Current simulation time in seconds type(RotParameterType), intent(in ) :: p !< AD parameters type(AD_ParameterType), intent(in ) :: p_AD !< AD parameters type(RotInputType), intent(in ) :: u !< AD Inputs at Time @@ -2648,12 +2660,17 @@ subroutine SetInputs(p, p_AD, u, m, indx, errStat, errMsg) ErrMsg = "" ! Disturbed inflow on blade (if tower shadow present) - call SetDisturbedInflow(p, p_AD, u, m, errStat, errMsg) + call SetDisturbedInflow(p, p_AD, u, m, errStat2, errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, RoutineName) + + if (p_AD%Wake_Mod /= WakeMod_FVW) then + + if (p_AD%SectAvg) then + call SetSectAvgInflow(t, p, p_AD, u, m, errStat2, errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, RoutineName) + endif - if (p_AD%WakeMod /= WakeMod_FVW) then ! This needs to extract the inputs from the AD data types (mesh) and massage them for the BEMT module - call SetInputsForBEMT(p, u, m, indx, errStat2, errMsg2) - call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) + call SetInputsForBEMT(p, p_AD, u, m, indx, errStat2, errMsg2) + call SetErrStat(errStat2, errMsg2, errStat, errMsg, RoutineName) endif end subroutine SetInputs @@ -2683,7 +2700,7 @@ subroutine SetDisturbedInflow(p, p_AD, u, m, errStat, errMsg) end do end if - if (p_AD%SkewMod == SkewMod_Orthogonal) then + if (p_AD%Skew_Mod == Skew_Mod_Orthogonal) then x_hat_disk = u%HubMotion%Orientation(1,:,1) do k=1,p%NumBlades @@ -2695,12 +2712,141 @@ subroutine SetDisturbedInflow(p, p_AD, u, m, errStat, errMsg) end subroutine SetDisturbedInflow +!> Sector Averaged (disturbed when tower influence on) inflow on the blade +!! Loop on blade nodes and computed a weighted sector average inflow at each node +subroutine SetSectAvgInflow(t, p, p_AD, u, m, errStat, errMsg) + real(DbKi), intent(in ) :: t !< Current simulation time in seconds + type(RotParameterType), intent(in ) :: p !< AD parameters + type(AD_ParameterType), intent(in ) :: p_AD !< AD parameters + type(RotInputType), intent(in ) :: u !< AD Inputs at Time + type(RotMiscVarType), intent(inout) :: m !< Misc/optimization variables + integer(IntKi), intent( out) :: errStat !< Error status of the operation + character(*), intent( out) :: errMsg !< Error message if ErrStat /= ErrID_None + ! local variables + real(R8Ki) :: R_li !< + real(ReKi) :: x_hat_disk(3) !< unit vector normal to disk along hub x axis + real(ReKi) :: r_A(3) !< Vector from global origin to blade node + real(ReKi) :: r_H(3) !< Vector from global origin to hub center + real(ReKi) :: r_S(3) !< Vector from global origin to point in sector + real(ReKi) :: rHS(3) !< Vector from rotor center to point in sector + real(ReKi) :: rHA(3) !< Vector from rotor center to blade node + real(ReKi) :: rHA_perp(3) !< Component of rHA perpendicular to x_hat_disk + real(ReKi) :: rHA_para(3) !< Component of rHA paralel to x_hat_disk + real(ReKi) :: rHA_perp_n !< Norm of rHA_perp + real(ReKi) :: e_r(3) !< Polar unit vector along rHA_perp + real(ReKi) :: e_t(3) !< Polar unit vector perpendicular to rHA_perp ("e_theta") + real(ReKi) :: temp_norm + real(ReKi) :: psi !< Azimuthal offset in the current sector, runs from -psi_bwd to psi_fwd + real(ReKi) :: dpsi !< Azimuthal increment + real(ReKi), allocatable :: SectPos(:,:)!< Points used to define a given sector (for a given blade node A) + real(ReKi), allocatable :: SectVel(:,:)!< Inflow velocity at a given sector (Undisturbed and then disturbed) + real(ReKi), allocatable :: SectAcc(:,:)!< Inflow velocity at a given sector (Undisturbed and then disturbed) + real(ReKi), allocatable :: SectWgt(:) !< Sector weights for velocity averaging + integer(intKi) :: j,k, ipsi + integer(intKi) :: errStat2 + character(ErrMsgLen) :: errMsg2 + character(*), parameter :: RoutineName = 'SetSectAvgInflow' + ! + errStat = ErrID_None + errMsg = "" + + if (.not. associated(p_AD%FlowField)) then + errStat2 = errID_Fatal + errMsg2 = 'FlowField should be allocated' + if (Failed()) return + endif + + ! Alloc and inits + call AllocAry(SectPos, 3, p_AD%SA_nPerSec, "SectPos", errStat2, errMsg2); if(Failed()) return + call AllocAry(SectVel, 3, p_AD%SA_nPerSec, "SectVel", errStat2, errMsg2); if(Failed()) return + call AllocAry(SectWgt, p_AD%SA_nPerSec, "SectWgt", errStat2, errMsg2); if(Failed()) return + if (allocated(SectAcc)) deallocate(SectAcc) ! IfW_FlowField_GetVelAcc some logic for Acc, so we ensure it's deallocated + SectVel = 0.0_ReKi + SectPos = 0.0_ReKi + if (p_AD%SA_Weighting == SA_Wgt_Uniform) then + SectWgt = 1.0_ReKi/p_AD%SA_nPerSec + else + errStat2 = errID_Fatal; errMsg2 = 'Sector averaging weighting (`SA_Weighting`) should be Uniform' + if (Failed()) return + endif + dpsi = (p_AD%SA_PsiFwd-p_AD%SA_PsiBwd)/(p_AD%SA_nPerSec-1) + + ! Hub + x_hat_disk = real(u%HubMotion%Orientation(1,:,1), ReKi) + r_H = u%HubMotion%Position(:,1) + u%HubMotion%TranslationDisp(:,1) + + ! --- Loop on blade nodes and computed a weighted sector average inflow at each node + do k=1,p%NumBlades + do j=1,p%NumBlNds + + ! --- Setup a polar coordinate system based on the current blade node + ! This is the same kind of calculations as the Calculate_MeshOrientation_Rel2Hub + r_A = u%BladeMotion(k)%Position(:,j) + u%BladeMotion(k)%TranslationDisp(:,j) + rHA = r_A - r_H + rHA_para = dot_product( x_hat_disk, rHA ) * x_hat_disk + rHA_perp = rHA - rHA_para + rHA_perp_n = TwoNorm( rHA_perp ) + + ! --- Create list of section points around the current blade node + if (EqualRealNos(rHA_perp_n, 0.0_ReKi)) then + ! We set all points to be the current one (likely the rotor center when no hub..) + do ipsi=1,p_AD%SA_nPerSec + SectPos(:, ipsi) = r_A + enddo + else + e_r = rHA_perp/rHA_perp_n ! Unit vector in "radial" coordinate + e_t = cross_product( x_hat_disk, e_r ) ! Unit vector in "tangential" coordinate + do ipsi=1,p_AD%SA_nPerSec + psi = p_AD%SA_PsiBwd + (ipsi-1)*dpsi + SectPos(:, ipsi) = (rHA_perp_n*cos(psi) * e_r + rHA_perp_n*sin(psi) * e_t) + rHA_para + r_H + enddo + endif + + ! --- Inflow on sector points + ! Undisturbed + call IfW_FlowField_GetVelAcc(p_AD%FlowField, 1, t, SectPos, SectVel, SectAcc, errStat=errStat2, errMsg=errMsg2); if(Failed()) return + ! --- Option 1 Disturbed inflow Before averaging - SectVel is modified in place + !if (p%TwrPotent /= TwrPotent_none .or. p%TwrShadow /= TwrShadow_none) then + ! call TwrInflArray(p, u, m, SectPos, SectVel, errStat2, errMsg2); if(Failed()) return + !endif + + ! --- Weighting and averaging + m%SectAvgInflow(1, j, k) = sum(SectVel(1,:)*SectWgt) + m%SectAvgInflow(2, j, k) = sum(SectVel(2,:)*SectWgt) + m%SectAvgInflow(3, j, k) = sum(SectVel(3,:)*SectWgt) + + ! --- Option 2 Disturbed after averaging + if (p%TwrPotent /= TwrPotent_none .or. p%TwrShadow /= TwrShadow_none) then + ! TODO use a "scalar" function or change the interface of TwrInfl. Waiting for Wind Inputs of AD to be removed from AD + call TwrInflArray( p, u, m, reshape(r_A, (/3,1/)), m%SectAvgInflow(:, j:j, k), errStat2, errMsg2); if(Failed()) return + endif + enddo + + enddo + + call CleanUp() +contains + subroutine CleanUp() + if(allocated(SectPos)) deallocate(SectPos) + if(allocated(SectVel)) deallocate(SectVel) + if(allocated(SectAcc)) deallocate(SectAcc) + if(allocated(SectWgt)) deallocate(SectWgt) + end subroutine + logical function Failed() + call SetErrStat(errStat2, errMsg2, errStat, errMsg, RoutineName) + Failed = errStat >= AbortErrLev + if (Failed) call CleanUp() + end function Failed +end subroutine SetSectAvgInflow + + !---------------------------------------------------------------------------------------------------------------------------------- !> This subroutine sets m%BEMT_u(indx). -subroutine SetInputsForBEMT(p, u, m, indx, errStat, errMsg) +subroutine SetInputsForBEMT(p, p_AD, u, m, indx, errStat, errMsg) type(RotParameterType), intent(in ) :: p !< AD parameters + type(AD_ParameterType), intent(in ) :: p_AD !< AD parameters type(RotInputType), intent(in ) :: u !< AD Inputs at Time type(RotMiscVarType), intent(inout) :: m !< Misc/optimization variables integer, intent(in ) :: indx !< index into m%BEMT_u array; must be 1 or 2 (but not checked here) @@ -2816,7 +2962,7 @@ subroutine SetInputsForBEMT(p, u, m, indx, errStat, errMsg) signofAngle = sign(1.0_ReKi,SkewVec(3)) endif - if (p%AeroBEM_Mod /= BEMMod_2D) then + if (p%BEM_Mod /= BEMMod_2D) then ! TODO m%BEMT_u(indx)%chi0 = sign( m%BEMT_u(indx)%chi0, signOfAngle ) endif end if @@ -2966,8 +3112,12 @@ subroutine SetInputsForBEMT(p, u, m, indx, errStat, errMsg) !.......................... do k=1,p%NumBlades do j=1,p%NumBlNds - ! Velocity in "l" or "w" system (depending) on AeroProjMod - tmp = m%DisturbedInflow(:,j,k) - u%BladeMotion(k)%TranslationVel(:,j) ! rel_V(j)_Blade(k) + if (p_AD%SectAvg) then + tmp = m%SectAvgInflow(:,j,k) - u%BladeMotion(k)%TranslationVel(:,j) ! rel_V(j)_Blade(k) + else + tmp = m%DisturbedInflow(:,j,k) - u%BladeMotion(k)%TranslationVel(:,j) ! rel_V(j)_Blade(k) + endif + ! Velocity in "p" or "w" system (depending) on AeroProjMod m%BEMT_u(indx)%Vx(j,k) = dot_product( tmp, m%orientationAnnulus(1,:,j,k) ) ! normal component (normal to the plane, not chord) of the inflow velocity of the jth node in the kth blade m%BEMT_u(indx)%Vy(j,k) = dot_product( tmp, m%orientationAnnulus(2,:,j,k) ) !+ TwoNorm(m%DisturbedInflow(:,j,k))*(sin()*sin(tilt)*)! tangential component (tangential to the plane, not chord) of the inflow velocity of the jth node in the kth blade m%BEMT_u(indx)%Vz(j,k) = dot_product( tmp, m%orientationAnnulus(3,:,j,k) ) ! radial component (tangential to the plane, not chord) of the inflow velocity of the jth node in the kth blade @@ -3304,10 +3454,22 @@ subroutine SetInputsForFVW(p, u, m, errStat, errMsg) ! Get disk average values and orientations ! NOTE: needed because it sets m%V_diskAvg and m%V_dot_x, needed by CalcOutput.. call DiskAvgValues(p%rotors(iR), u(tIndx)%rotors(iR), m%rotors(iR), x_hat_disk) ! also sets m%V_diskAvg and m%V_dot_x + + ! Compute Orientation similar to BEM, only to have consistent outputs... + ! TODO TODO TODO All this below is mostly a calcOutput thing, we should move it somewhere else! + ! orientation annulus is only used for Outputs with OLAF, same for pitch and azimuth if (p%rotors(iR)%AeroProjMod==APM_BEM_NoSweepPitchTwist) then - call Calculate_MeshOrientation_NoSweepPitchTwist(p%rotors(iR),u(tIndx)%rotors(iR), m%rotors(iR), thetaBladeNds,ErrStat=ErrStat2,ErrMsg=ErrMsg2) ! sets m%orientationAnnulus, m%Curve + call Calculate_MeshOrientation_NoSweepPitchTwist(p%rotors(iR), u(tIndx)%rotors(iR), m%rotors(iR), thetaBladeNds,ErrStat=ErrStat2,ErrMsg=ErrMsg2) ! sets m%orientationAnnulus, m%Curve + + elseif (p%rotors(iR)%AeroProjMod==APM_BEM_Polar) then + do k=1,p%rotors(iR)%numBlades + call Calculate_MeshOrientation_Rel2Hub(u(tIndx)%rotors(iR)%BladeMotion(k), u(tIndx)%rotors(iR)%HubMotion, x_hat_disk, m%rotors(iR)%orientationAnnulus(:,:,:,k)) + enddo + else if (p%rotors(iR)%AeroProjMod==APM_LiftingLine) then call Calculate_MeshOrientation_LiftingLine (p%rotors(iR),u(tIndx)%rotors(iR), m%rotors(iR), thetaBladeNds,ErrStat=ErrStat2,ErrMsg=ErrMsg2) ! sets m%orientationAnnulus, m%Curve + else + call SetErrStat(ErrID_Fatal, 'Aero Projection Method not implemented' ,ErrStat, ErrMsg, RoutineName) endif call StorePitchAndAzimuth(p%rotors(iR), u(tIndx)%rotors(iR), m%rotors(iR), ErrStat2, ErrMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) @@ -3416,7 +3578,7 @@ subroutine SetOutputsFromBEMT( p, u, m, y ) real(reki) :: c ! local chord length real(reki) :: aoa ! local angle of attack real(reki) :: Cl,Cd,Cm ! local airfoil lift, drag and pitching moment coefficients - real(reki) :: Cn,Ct ! local airfoil normal and tangential force coefficients + real(reki) :: Cxa,Cya ! local airfoil normal and tangential force coefficients do k=1,p%NumBlades @@ -3427,14 +3589,14 @@ subroutine SetOutputsFromBEMT( p, u, m, y ) Cl = m%BEMT_y%cl(j,k) Cd = m%BEMT_y%cd(j,k) Cm = m%BEMT_y%cm(j,k) - Cn = Cl*cos(aoa) + Cd*sin(aoa) - Ct = -Cl*sin(aoa) + Cd*cos(aoa) ! NOTE: this is not Ct but Cy_a (y_a going towards the TE) + Cxa = Cl*cos(aoa) + Cd*sin(aoa) + Cya = -Cl*sin(aoa) + Cd*cos(aoa) ! Dimensionalize the aero forces and moment q = 0.5 * p%airDens * m%BEMT_y%Vrel(j,k)**2 ! dynamic pressure of the jth node in the kth blade c = p%BEMT%chord(j,k) - forceAirfoil(1) = Cn * q * c - forceAirfoil(2) = Ct * q * c + forceAirfoil(1) = Cxa * q * c + forceAirfoil(2) = Cya * q * c forceAirfoil(3) = 0.0_reki momentAirfoil(1) = 0.0_reki momentAirfoil(2) = 0.0_reki @@ -3667,15 +3829,11 @@ SUBROUTINE ValidateInputData( InitInp, InputFileData, NumBl, ErrStat, ErrMsg ) ! end do if (InputFileData%DTAero <= 0.0) call SetErrStat ( ErrID_Fatal, 'DTAero must be greater than zero.', ErrStat, ErrMsg, RoutineName ) - if (InputFileData%WakeMod /= WakeMod_None .and. InputFileData%WakeMod /= WakeMod_BEMT .and. InputFileData%WakeMod /= WakeMod_DBEMT .and. InputFileData%WakeMod /= WakeMod_FVW) then + if (InputFileData%Wake_Mod /= WakeMod_None .and. InputFileData%Wake_Mod /= WakeMod_BEMT .and. InputFileData%Wake_Mod /= WakeMod_FVW) then call SetErrStat ( ErrID_Fatal, 'WakeMod must be '//trim(num2lstr(WakeMod_None))//' (none), '//trim(num2lstr(WakeMod_BEMT))//' (BEMT), '// & - trim(num2lstr(WakeMod_DBEMT))//' (DBEMT), or '//trim(num2lstr(WakeMod_FVW))//' (FVW).',ErrStat, ErrMsg, RoutineName ) + ' or '//trim(num2lstr(WakeMod_FVW))//' (FVW).',ErrStat, ErrMsg, RoutineName ) end if - if (InputFileData%AFAeroMod /= AFAeroMod_Steady .and. InputFileData%AFAeroMod /= AFAeroMod_BL_unsteady) then - call SetErrStat ( ErrID_Fatal, 'AFAeroMod must be '//trim(num2lstr(AFAeroMod_Steady))//' (steady) or '//& - trim(num2lstr(AFAeroMod_BL_unsteady))//' (Beddoes-Leishman unsteady).', ErrStat, ErrMsg, RoutineName ) - end if if (InputFileData%TwrPotent /= TwrPotent_none .and. InputFileData%TwrPotent /= TwrPotent_baseline .and. InputFileData%TwrPotent /= TwrPotent_Bak) then call SetErrStat ( ErrID_Fatal, 'TwrPotent must be 0 (none), 1 (baseline potential flow), or 2 (potential flow with Bak correction).', ErrStat, ErrMsg, RoutineName ) end if @@ -3704,6 +3862,7 @@ SUBROUTINE ValidateInputData( InitInp, InputFileData, NumBl, ErrStat, ErrMsg ) if ( maxval(InputFileData%rotors(iR)%TwrTI) > 0.4 .and. maxval(InputFileData%rotors(iR)%TwrTI) < 1.0) call SetErrStat ( ErrID_Warn, 'The turbulence intensity for the Eames tower shadow model above 0.4 may return unphysical results. Interpret with caution.', ErrStat, ErrMsg, RoutineName ) enddo endif + if (Failed()) return if (InitInp%MHK == MHK_None .and. InputFileData%CavitCheck) call SetErrStat ( ErrID_Fatal, 'A cavitation check can only be performed for an MHK turbine.', ErrStat, ErrMsg, RoutineName ) if (InitInp%MHK == MHK_None .and. InputFileData%Buoyancy) call SetErrStat ( ErrID_Fatal, 'Buoyancy can only be calculated for an MHK turbine.', ErrStat, ErrMsg, RoutineName ) @@ -3720,21 +3879,34 @@ SUBROUTINE ValidateInputData( InitInp, InputFileData, NumBl, ErrStat, ErrMsg ) - ! BEMT/DBEMT inputs - ! bjj: these checks should probably go into BEMT where they are used... - if (InputFileData%WakeMod /= WakeMod_none .and. InputFileData%WakeMod /= WakeMod_FVW) then + ! NOTE: this check is done here because it is used for all kind of Wake Mod + if (.not.any(InputFileData%BEM_Mod == (/BEMMod_2D, BEMMod_3D/))) call Fatal('BEM_Mod must be 1 or 2.') + + ! --- BEMT/DBEMT inputs + ! bjj: these checks should probably go into BEMT where they are used... + if (InputFileData%Wake_Mod == WakeMod_BEMT) then if ( InputFileData%MaxIter < 1 ) call SetErrStat( ErrID_Fatal, 'MaxIter must be greater than 0.', ErrStat, ErrMsg, RoutineName ) if ( InputFileData%IndToler < 0.0 .or. EqualRealNos(InputFileData%IndToler, 0.0_ReKi) ) & call SetErrStat( ErrID_Fatal, 'IndToler must be greater than 0.', ErrStat, ErrMsg, RoutineName ) - if ( InputFileData%SkewMod /= SkewMod_Orthogonal .and. InputFileData%SkewMod /= SkewMod_Uncoupled .and. InputFileData%SkewMod /= SkewMod_PittPeters) & ! .and. InputFileData%SkewMod /= SkewMod_Coupled ) - call SetErrStat( ErrID_Fatal, 'SkewMod must be 1, or 2. Option 3 will be implemented in a future version.', ErrStat, ErrMsg, RoutineName ) + if (.not.any(InputFileData%Skew_Mod == (/Skew_Mod_Orthogonal, Skew_Mod_None, Skew_Mod_Active/))) call Fatal('Skew_Mod must be -1, 0, or 1.') + if (.not.any(InputFileData%SkewRedistr_Mod == (/SkewRedistrMod_None, SkewRedistrMod_PittPeters/))) call Fatal('SkewRedistr_Mod should be 0 or 1') + + if ( InputFileData%SectAvg) then + if (InputFileData%SA_Weighting /= SA_Wgt_Uniform) call Fatal('SectAvgWeighting should be Uniform (=1) for now.') + if (InputFileData%SA_nPerSec <= 1) call Fatal('SectAvgNPoints must be >=1') + if (InputFileData%SA_PsiBwd > 0) call Fatal('SectAvgPsiBwd must be negative') + if (InputFileData%SA_PsiFwd < 0) call Fatal('SectAvgPsiFwd must be positive') + if (InputFileData%SA_PsiFwd <= InputFileData%SA_PsiBwd ) call Fatal('SectAvgPsiFwd must be strictly higher than SA_PsiBwd') + endif + ! Good to return once in a while.. + if (Failed()) return end if !BEMT/DBEMT checks - if ( InputFileData%CavitCheck .and. InputFileData%AFAeroMod == AFAeroMod_BL_unsteady) then + if ( InputFileData%CavitCheck .and. InputFileData%UA_Mod >0) then call SetErrStat( ErrID_Fatal, 'Cannot use unsteady aerodynamics module with a cavitation check', ErrStat, ErrMsg, RoutineName ) end if @@ -3799,6 +3971,7 @@ SUBROUTINE ValidateInputData( InitInp, InputFileData, NumBl, ErrStat, ErrMsg ) end do ! k=blades end if end do ! iR rotor + if (Failed()) return ! ............................. ! check tower mesh data: @@ -3841,6 +4014,7 @@ SUBROUTINE ValidateInputData( InitInp, InputFileData, NumBl, ErrStat, ErrMsg ) end if end do ! iR rotor + if (Failed()) return @@ -3892,6 +4066,7 @@ SUBROUTINE ValidateInputData( InitInp, InputFileData, NumBl, ErrStat, ErrMsg ) enddo ! iR end if + if (Failed()) return if ( ( InputFileData%NBlOuts < 0_IntKi ) .OR. ( InputFileData%NBlOuts > 9_IntKi ) ) then @@ -3911,6 +4086,7 @@ SUBROUTINE ValidateInputData( InitInp, InputFileData, NumBl, ErrStat, ErrMsg ) end do ! iR, rotor end if + if (Failed()) return !.................. ! Tail fin checks @@ -3925,23 +4101,23 @@ SUBROUTINE ValidateInputData( InitInp, InputFileData, NumBl, ErrStat, ErrMsg ) endif endif enddo ! iR, rotor + if (Failed()) return !.................. ! check for linearization !.................. if (InitInp%Linearize) then - if (InputFileData%AFAeroMod /= AFAeroMod_Steady) then - if (InputFileData%UAMod /= UA_HGM .and. InputFileData%UAMod /= UA_HGMV .and. InputFileData%UAMod /= UA_OYE) then - call SetErrStat( ErrID_Fatal, 'When AFAeroMod=2, UAMod must be 4, 5, or 6 for linearization. Set AFAeroMod=1, or, set UAMod=4, 5, or 6.', ErrStat, ErrMsg, RoutineName ) - end if + + if (InputFileData%Wake_Mod /= WakeMod_None .and. InputFileData%Wake_Mod /= WakeMod_BEMT) then + call SetErrStat( ErrID_Fatal, 'WakeMod must be 0 or 1 for linearization.', ErrStat, ErrMsg, RoutineName ) + endif + + if (InputFileData%UA_Mod /= UA_None .and. InputFileData%UA_Mod /= UA_HGM .and. InputFileData%UA_Mod /= UA_HGMV .and. InputFileData%UA_Mod /= UA_OYE) then + call SetErrStat( ErrID_Fatal, 'UA_Mod must be 0, 4, 5, or 6 for linearization.', ErrStat, ErrMsg, RoutineName ) end if - - if (InputFileData%WakeMod == WakeMod_FVW) then !bjj: note: among other things, WriteOutput values will not be calculated properly in AD Jacobians if FVW this is allowed - call SetErrStat( ErrID_Fatal, 'FVW cannot currently be used for linearization. Set WakeMod=0 or WakeMod=1.', ErrStat, ErrMsg, RoutineName ) - else if (InputFileData%WakeMod == WakeMod_DBEMT) then - if (InputFileData%DBEMT_Mod /= DBEMT_cont_tauConst) then - call SetErrStat( ErrID_Fatal, 'DBEMT requires the continuous formulation with constant tau1 for linearization. Set DBEMT_Mod=3 or set WakeMod to 0 or 1.', ErrStat, ErrMsg, RoutineName ) - end if + + if (InputFileData%DBEMT_Mod /= DBEMT_None .and. InputFileData%DBEMT_Mod /= DBEMT_cont_tauConst) then + call SetErrStat( ErrID_Fatal, 'DBEMT Mod must be 0 or 3 (continuous formulation with constant tau1) for linearization. Set DBEMT_Mod=0,3.', ErrStat, ErrMsg, RoutineName ) end if end if @@ -3951,6 +4127,10 @@ SUBROUTINE Fatal(ErrMsg_in) character(*), intent(in) :: ErrMsg_in call SetErrStat(ErrID_Fatal, ErrMsg_in, ErrStat, ErrMsg, RoutineName) END SUBROUTINE Fatal + + logical function Failed() + Failed = ErrStat >= AbortErrLev + end function Failed END SUBROUTINE ValidateInputData !---------------------------------------------------------------------------------------------------------------------------------- @@ -3994,7 +4174,7 @@ SUBROUTINE Init_AFIparams( InputFileData, p_AFI, UnEc, ErrStat, ErrMsg ) IF (.not. InputFileData%UseBlCm) AFI_InitInputs%InCol_Cm = 0 ! Don't try to use Cm if flag set to false AFI_InitInputs%InCol_Cpmin = InputFileData%InCol_Cpmin AFI_InitInputs%AFTabMod = InputFileData%AFTabMod !AFITable_1 - AFI_InitInputs%UA_f_cn = (InputFileData%UAMod /= UA_HGM).and.(InputFileData%UAMod /= UA_OYE) ! HGM and OYE use the separation function based on cl instead of cn + AFI_InitInputs%UA_f_cn = (InputFileData%UA_Mod /= UA_HGM).and.(InputFileData%UA_Mod /= UA_OYE) ! HGM and OYE use the separation function based on cl instead of cn ! Call AFI_Init to read in and process the airfoil files. ! This includes creating the spline coefficients to be used for interpolation. @@ -4155,6 +4335,7 @@ SUBROUTINE Init_BEMTmodule( InputFileData, RotInputFileData, u_AD, u, p, p_AD, x integer(IntKi) :: ErrStat2 character(ErrMsgLen) :: ErrMsg2 character(*), parameter :: RoutineName = 'Init_BEMTmodule' + character(1024) :: Label ! note here that each blade is required to have the same number of nodes @@ -4168,12 +4349,21 @@ SUBROUTINE Init_BEMTmodule( InputFileData, RotInputFileData, u_AD, u, p, p_AD, x InitInp%airDens = InputFileData%AirDens InitInp%kinVisc = InputFileData%KinVisc - InitInp%skewWakeMod = InputFileData%SkewMod + ! --- Skew + InitInp%skewWakeMod = InputFileData%Skew_Mod + InitInp%skewRedistrMod = InputFileData%SkewRedistr_Mod InitInp%yawCorrFactor = InputFileData%SkewModFactor + InitInp%MomentumCorr = InputFileData%SkewMomCorr + ! Safety + if (InputFileData%Skew_Mod /= Skew_Mod_Active) then + InitInp%skewRedistrMod = SkewRedistrMod_None + InitInp%MomentumCorr = .False. + endif + ! --- Algo InitInp%aTol = InputFileData%IndToler InitInp%useTipLoss = InputFileData%TipLoss InitInp%useHubLoss = InputFileData%HubLoss - InitInp%useInduction = InputFileData%WakeMod /= WakeMod_none + InitInp%useInduction = InputFileData%Wake_Mod == WakeMod_BEMT InitInp%useTanInd = InputFileData%TanInd InitInp%useAIDrag = InputFileData%AIDrag InitInp%useTIDrag = InputFileData%TIDrag @@ -4274,39 +4464,52 @@ SUBROUTINE Init_BEMTmodule( InputFileData, RotInputFileData, u_AD, u, p, p_AD, x end do InitInp%UA_Flag = p_AD%UA_Flag - InitInp%UAMod = InputFileData%UAMod + InitInp%UAMod = InputFileData%UA_Mod InitInp%Flookup = InputFileData%Flookup InitInp%a_s = InputFileData%SpdSound - InitInp%MomentumCorr = .FALSE. ! TODO EB InitInp%SumPrint = InputFileData%SumPrint InitInp%RootName = p%RootName - InitInp%BEM_Mod = p%AeroBEM_Mod - + InitInp%BEM_Mod = InputFileData%BEM_Mod + p%BEM_Mod = InputFileData%BEM_Mod ! TODO try to get rid of me - if (p%AeroBEM_Mod==-1) then - !call WrSCr('WARNING: AeroDyn: BEM_Mod is -1, using default BEM_Mod based on projection') - if (p%AeroProjMod == APM_BEM_NoSweepPitchTwist) then - InitInp%BEM_Mod = BEMMod_2D - else if (p%AeroProjMod == APM_BEM_Polar) then - InitInp%BEM_Mod = BEMMod_3D - else - InitInp%BEM_Mod = -1 - call SetErrStat(ErrID_Fatal, "AeroProjMod needs to be 1 or 2 when used with BEM", ErrStat, ErrMsg, RoutineName) - endif + ! --- Print BEM formulation to screen + Label = '' + if (p%AeroProjMod==APM_BEM_NoSweepPitchTwist) then + Label='Projection: legacy (NoSweepPitchTwist)' + elseif (p%AeroProjMod==APM_BEM_Polar) then + Label='Projection: Polar' + elseif (p%AeroProjMod==APM_LiftingLine) then + Label='Projection: Lifting Line' + else + ! Normally we wouldn't want to do a print or STOP, but we + ! should never get here unless a programmer made a mistake. + ! I'll leave this as is for now. - ADP + print*,'Invalid projection method' + STOP endif - p%AeroBEM_Mod = InitInp%BEM_Mod ! Very important, for consistency - !call WrScr(' AeroDyn: projMod: '//trim(num2lstr(p%AeroProjMod))//', BEM_Mod:'//trim(num2lstr(InitInp%BEM_Mod))) + if (InitInp%BEM_Mod==BEMMod_2D) then + Label = trim(Label)//', BEM: legacy (2D)' + elseif (InitInp%BEM_Mod==BEMMod_3D) then + Label = trim(Label)//', BEM: polar (3D)' + else + print*,'Invalid BEM method' + STOP + endif + if (InitInp%MomentumCorr) then + Label = trim(Label)//', MomentumCorrection' + endif + if (p_AD%SectAvg) then + Label = trim(Label)//', Sector Average' + endif + call WrScr(' '//trim(Label)) + ! remove the ".AD" from the RootName k = len_trim(InitInp%RootName) if (k>3) then InitInp%RootName = InitInp%RootName(1:k-3) end if - if (InputFileData%WakeMod == WakeMod_DBEMT) then - InitInp%DBEMT_Mod = InputFileData%DBEMT_Mod - else - InitInp%DBEMT_Mod = DBEMT_none - end if + InitInp%DBEMT_Mod = InputFileData%DBEMT_Mod InitInp%tau1_const = InputFileData%tau1_const if (ErrStat >= AbortErrLev) then @@ -4458,7 +4661,7 @@ SUBROUTINE Init_OLAF( InputFileData, u_AD, u, p, x, xd, z, OtherState, m, ErrSta ! Unsteady Aero Data InitInp%UA_Flag = p%UA_Flag - InitInp%UAMod = InputFileData%UAMod + InitInp%UAMod = InputFileData%UA_Mod InitInp%Flookup = InputFileData%Flookup InitInp%a_s = InputFileData%SpdSound InitInp%SumPrint = InputFileData%SumPrint @@ -4535,7 +4738,7 @@ SUBROUTINE TFin_CalcOutput(p, p_AD, u, m, y, ErrStat, ErrMsg ) elseif(p%TFin%TFinIndMod==TFinIndMod_rotavg) then ! TODO TODO - print*,'TODO TailFin: compute rotor average induced velocity' + call WrScr('TODO TailFin: compute rotor average induced velocity') V_ind = 0.0_ReKi else @@ -5329,8 +5532,8 @@ SUBROUTINE Rot_JacobianPInput( t, u, p, p_AD, x, xd, z, OtherState, y, m, m_AD, ! get OP values here (i.e., set inputs for BEMT): - if ( p%FrozenWake ) then - call SetInputs(p, p_AD, u, m, indx, errStat2, errMsg2) + if ( p%DBEMT_Mod == DBEMT_frozen ) then + call SetInputs(t, p, p_AD, u, m, indx, errStat2, errMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) ! we shouldn't have any errors about allocating memory here so I'm not going to return-on-error until later ! compare m%BEMT_y arguments with call to BEMT_CalcOutput @@ -5351,7 +5554,7 @@ SUBROUTINE Rot_JacobianPInput( t, u, p, p_AD, x, xd, z, OtherState, y, m, m_AD, ! initialize x_init so that we get accurrate values for first step if (.not. OtherState%BEMT%nodesInitialized ) then - call SetInputs(p, p_AD, u, m, indx, errStat2, errMsg2) + call SetInputs(t, p, p_AD, u, m, indx, errStat2, errMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) call BEMT_InitStates(t, m%BEMT_u(indx), p%BEMT, x_init%BEMT, xd%BEMT, z%BEMT, OtherState_init%BEMT, m%BEMT, p_AD%AFI, ErrStat2, ErrMsg2 ) ! changes values only if states haven't been initialized @@ -5414,7 +5617,7 @@ SUBROUTINE Rot_JacobianPInput( t, u, p, p_AD, x, xd, z, OtherState, y, m, m_AD, !call AD_UpdateStates( t, 1, (/u_perturb/), (/t/), p, x_copy, xd_copy, z_copy, OtherState_copy, m, errStat2, errMsg2 ) ! call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) !bjj: this is what we want to do instead of the overkill of calling AD_UpdateStates - call SetInputs(p, p_AD, u_perturb, m, indx, errStat2, errMsg2) + call SetInputs(t, p, p_AD, u_perturb, m, indx, errStat2, errMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) ! we shouldn't have any errors about allocating memory here so I'm not going to return-on-error until later call UpdatePhi( m%BEMT_u(indx), p%BEMT, z_copy%BEMT%phi, p_AD%AFI, m%BEMT, OtherState_copy%BEMT%ValidPhi, errStat2, errMsg2 ) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) ! we shouldn't have any errors about allocating memory here so I'm not going to return-on-error until later @@ -5437,7 +5640,7 @@ SUBROUTINE Rot_JacobianPInput( t, u, p, p_AD, x, xd, z, OtherState, y, m, m_AD, ! get updated z%phi values: !call RotUpdateStates( t, 1, (/u_perturb/), (/t/), p, x_copy, xd_copy, z_copy, OtherState_copy, m, errStat2, errMsg2 ) ! call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) - call SetInputs(p, p_AD, u_perturb, m, indx, errStat2, errMsg2) + call SetInputs(t, p, p_AD, u_perturb, m, indx, errStat2, errMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) ! we shouldn't have any errors about allocating memory here so I'm not going to return-on-error until later call UpdatePhi( m%BEMT_u(indx), p%BEMT, z_copy%BEMT%phi, p_AD%AFI, m%BEMT, OtherState_copy%BEMT%ValidPhi, errStat2, errMsg2 ) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) ! we shouldn't have any errors about allocating memory here so I'm not going to return-on-error until later @@ -5647,8 +5850,8 @@ SUBROUTINE RotJacobianPContState( t, u, p, p_AD, x, xd, z, OtherState, y, m, m_A ErrMsg = '' - if ( p%FrozenWake ) then - call SetInputs(p, p_AD, u, m, indx, errStat2, errMsg2) + if ( p%DBEMT_Mod == DBEMT_frozen ) then + call SetInputs(t, p, p_AD, u, m, indx, errStat2, errMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) ! compare arguments with call to BEMT_CalcOutput @@ -5672,7 +5875,7 @@ SUBROUTINE RotJacobianPContState( t, u, p, p_AD, x, xd, z, OtherState, y, m, m_A ! initialize x_init so that we get accurrate values for if (.not. OtherState%BEMT%nodesInitialized ) then - call SetInputs(p, p_AD, u, m, indx, errStat2, errMsg2) + call SetInputs(t, p, p_AD, u, m, indx, errStat2, errMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) call BEMT_InitStates(t, m%BEMT_u(indx), p%BEMT, x_init%BEMT, xd%BEMT, z%BEMT, OtherState_init%BEMT, m%BEMT, p_AD%AFI, ErrStat2, ErrMsg2 ) ! changes values only if states haven't been initialized @@ -6019,13 +6222,13 @@ SUBROUTINE RotJacobianPConstrState( t, u, p, p_AD, x, xd, z, OtherState, y, m, m ! get OP values here: !call AD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, ErrStat2, ErrMsg2 ) ! (bjj: is this necessary? if not, still need to get BEMT inputs) - call SetInputs(p, p_AD, u, m, indx, errStat2, errMsg2) + call SetInputs(t, p, p_AD, u, m, indx, errStat2, errMsg2) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) ! we shouldn't have any errors about allocating memory here so I'm not going to return-on-error until later call BEMT_CopyInput( m%BEMT_u(indx), m%BEMT_u(op_indx), MESH_UPDATECOPY, ErrStat2, ErrMsg2) ! copy the BEMT OP inputs to a temporary location that won't be overwritten call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) ! we shouldn't have any errors about allocating memory here so I'm not going to return-on-error until later - if ( p%FrozenWake ) then + if ( p%DBEMT_Mod == DBEMT_frozen ) then ! compare arguments with call to BEMT_CalcOutput call computeFrozenWake(m%BEMT_u(op_indx), p%BEMT, m%BEMT_y, m%BEMT ) m%BEMT%UseFrozenWake = .true. @@ -7158,7 +7361,7 @@ SUBROUTINE Init_Jacobian( InputFileData, p, p_AD, u, y, m, InitOut, ErrStat, Err call Init_Jacobian_y( p, p_AD, y, InitOut, ErrStat, ErrMsg) ! these matrices will be needed for linearization with frozen wake feature - if (p%FrozenWake) then + if ( p%DBEMT_Mod == DBEMT_frozen ) then call AllocAry(m%BEMT%AxInd_op,p%NumBlNds,p%numBlades,'m%BEMT%AxInd_op', ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) call AllocAry(m%BEMT%TnInd_op,p%NumBlNds,p%numBlades,'m%BEMT%TnInd_op', ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) end if diff --git a/modules/aerodyn/src/AeroDyn_AllBldNdOuts_IO.f90 b/modules/aerodyn/src/AeroDyn_AllBldNdOuts_IO.f90 index 2bc7877acf..27a3e3899b 100644 --- a/modules/aerodyn/src/AeroDyn_AllBldNdOuts_IO.f90 +++ b/modules/aerodyn/src/AeroDyn_AllBldNdOuts_IO.f90 @@ -301,7 +301,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx nNd = p%NumBlNds R_li => m%R_li ! inertial to local-polar R_wi => m%orientationAnnulus ! inertial to without-sweep-pitch-twist or orientation annulus (TODO: deprecate me) - if (p_AD%WakeMod == WakeMod_FVW) W2B => p_AD%FVW%Bld2Wings(iRot, :) ! From Wing index to blade index + if (p_AD%Wake_Mod == WakeMod_FVW) W2B => p_AD%FVW%Bld2Wings(iRot, :) ! From Wing index to blade index ! Initialize some things ErrMsg = '' @@ -406,7 +406,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx ! Axial and tangential induced wind velocity ! TODO use m%Vind_i and R_wi CASE ( BldNd_Vindx ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then do iB=1,nB do iNd=1,nNd y%WriteOutput(iOut) = - m%BEMT_u(Indx)%Vx(iNd,iB) * m%BEMT_y%axInduction( iNd,iB) @@ -424,7 +424,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx endif CASE ( BldNd_Vindy ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = m%BEMT_u(Indx)%Vy(iNd,iB) * m%BEMT_y%tanInduction(iNd,iB) @@ -461,7 +461,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx ! TODO: Vrel, DynP, Re, Ma - should be unified across lifting-line implementations. Vrel should be computed based on velocities in (a)-system ! Relative wind speed CASE ( BldNd_VRel ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = m%BEMT_y%Vrel(iNd,iB) @@ -480,7 +480,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx ! Dynamic pressure CASE ( BldNd_DynP ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = 0.5 * p%airDens * m%BEMT_y%Vrel(iNd,iB)**2 @@ -499,7 +499,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx ! Reynolds number (in millions) CASE ( BldNd_Re ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = p%BEMT%chord(iNd,iB) * m%BEMT_y%Vrel(iNd,iB) / p%KinVisc / 1.0E6 @@ -518,7 +518,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx ! Mach number CASE ( BldNd_M ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = m%BEMT_y%Vrel(iNd,iB) / p%SpdSound @@ -538,7 +538,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx ! Axial and tangential induction factors CASE ( BldNd_AxInd ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = m%BEMT_y%axInduction(iNd,iB) @@ -556,7 +556,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx endif CASE ( BldNd_TnInd ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = m%BEMT_y%tanInduction(iNd,iB) @@ -575,7 +575,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx ! Quasi-steady Axial and tangential induction factors CASE ( BldNd_AxInd_qs ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = m%BEMT_y%axInduction_qs(iNd,iB) @@ -593,7 +593,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx endif CASE ( BldNd_TnInd_qs ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = m%BEMT_y%tanInduction_qs(iNd,iB) @@ -610,10 +610,10 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx END DO endif - + ! AoA, pitch+twist angle, inflow angle, and curvature angle CASE ( BldNd_Alpha ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd ! TODO Change this @@ -632,7 +632,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx endif CASE ( BldNd_Theta ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = m%BEMT_u(Indx)%theta(iNd,iB)*R2D @@ -650,7 +650,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx endif CASE ( BldNd_Phi ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = m%BEMT_y%phi(iNd,iB)*R2D @@ -668,7 +668,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx endif CASE ( BldNd_Curve ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = m%Curve(iNd,iB)*R2D @@ -687,7 +687,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx endif CASE ( BldNd_Toe ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = m%BEMT_u(Indx)%toeAngle(iNd,iB)*R2D @@ -708,7 +708,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx ! Unsteady lift force, drag force, pitching moment coefficients ! TODO this should be somehow unified across lifting-line implementations CASE ( BldNd_Cl ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = m%BEMT_y%Cl(iNd,iB) @@ -726,7 +726,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx endif CASE ( BldNd_Cd ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = m%BEMT_y%Cd(iNd,iB) @@ -744,7 +744,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx endif CASE ( BldNd_Cm ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = m%BEMT_y%Cm(iNd,iB) @@ -764,7 +764,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx ! Normal force (to plane), tangential force (to plane) coefficients ! TODO deprecate CASE ( BldNd_Cx ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = m%BEMT_y%Cx(iNd,iB) @@ -782,7 +782,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx endif CASE ( BldNd_Cy ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = m%BEMT_y%Cy(iNd,iB) @@ -801,7 +801,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx ! Normal force (to chord), and tangential force (to chord) coefficients CASE ( BldNd_Cn ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd ct=cos(m%BEMT_u(Indx)%theta(iNd,iB)) @@ -823,7 +823,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx endif CASE ( BldNd_Ct ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd ct=cos(m%BEMT_u(Indx)%theta(iNd,iB)) @@ -847,7 +847,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx ! Lift force, drag force, pitching moment CASE ( BldNd_Fl ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd cp=cos(m%BEMT_y%phi(iNd,iB)) @@ -869,7 +869,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx endif CASE ( BldNd_Fd ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd cp=cos(m%BEMT_y%phi(iNd,iB)) @@ -900,7 +900,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx ! Normal force (to chord), and tangential force (to chord) per unit length !CASE( BldNd_Fn ); do iB=1,nB; do iNd=1,nNd; y%WriteOutput(iOut) = dot_product( y%BladeLoad(iB)%Force (:, iNd), u%BladeMotion(iB)%Orientation(1,:,iNd)); iOut = iOut + 1; enddo;enddo CASE ( BldNd_Fn ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd ct=cos(m%BEMT_u(Indx)%theta(iNd,iB)) @@ -923,7 +923,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx !CASE( BldNd_Ft ); do iB=1,nB; do iNd=1,nNd; y%WriteOutput(iOut) = -dot_product( y%BladeLoad(iB)%Force (:, iNd), u%BladeMotion(iB)%Orientation(2,:,iNd)); iOut = iOut + 1; enddo;enddo CASE ( BldNd_Ft ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd ct=cos(m%BEMT_u(Indx)%theta(iNd,iB)) @@ -992,7 +992,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx ! TODO: remove me, Vx, Vy can be computed from other outputs (and they are in legacy coordinate system) CASE ( BldNd_Vx ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = m%BEMT_u(Indx)%Vx(iNd,iB) @@ -1010,7 +1010,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx endif CASE ( BldNd_Vy ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = m%BEMT_u(Indx)%Vy(iNd,iB) @@ -1028,7 +1028,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx endif CASE ( BldNd_GeomPhi ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then if (allocated(OtherState%BEMT%ValidPhi)) then DO iB=1,nB DO iNd=1,nNd @@ -1059,7 +1059,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx endif CASE ( BldNd_chi ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,nNd y%WriteOutput(iOut) = m%BEMT_y%chi(iNd,iB)*R2D @@ -1078,7 +1078,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx CASE ( BldNd_UA_Flag ) IF (p_AD%UA_Flag) THEN - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,u%BladeMotion(iB)%NNodes y%WriteOutput(iOut) = m%BEMT%UA%weight(iNd, iB) @@ -1119,7 +1119,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx compIndx = 5 END SELECT - !if (p_AD%WakeMod /= WakeMod_FVW) then + !if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,u%BladeMotion(iB)%NNodes y%WriteOutput(iOut) = x%BEMT%UA%element(iNd, iB)%x(compIndx) @@ -1148,7 +1148,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx ! CpMin CASE ( BldNd_CpMin ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,u%BladeMotion(iB)%NNodes y%WriteOutput(iOut) = m%BEMT_y%Cpmin(iNd,iB) @@ -1171,7 +1171,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx ! circulation on blade CASE ( BldNd_Gam ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,u%BladeMotion(iB)%NNodes y%WriteOutput(iOut) = 0.5_ReKi * p%BEMT%chord(iNd,iB) * m%BEMT_y%Vrel(iNd,iB) * m%BEMT_y%Cl(iNd,iB) ! "Gam" [m^2/s] @@ -1194,7 +1194,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx ! TODO this should be provided by all lifting-line codes ! Cl_Static CASE ( BldNd_Cl_qs ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,u%BladeMotion(iB)%NNodes !NOT available in BEMT/DBEMT yet @@ -1214,7 +1214,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx ! Cd_Static CASE ( BldNd_Cd_qs ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,u%BladeMotion(iB)%NNodes !NOT available in BEMT/DBEMT yet @@ -1234,7 +1234,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, Indx ! Cm_Static CASE ( BldNd_Cm_qs ) - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then DO iB=1,nB DO iNd=1,u%BladeMotion(iB)%NNodes !NOT available in BEMT/DBEMT yet @@ -1531,7 +1531,7 @@ SUBROUTINE BldNdOuts_SetOutParam(BldNd_OutList, p, p_AD, ErrStat, ErrMsg ) end if ! The following are valid only for BEMT/DBEMT - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then InvalidOutput( BldNd_Cl_qs ) = .true. InvalidOutput( BldNd_Cd_qs ) = .true. InvalidOutput( BldNd_Cm_qs ) = .true. @@ -1544,7 +1544,7 @@ SUBROUTINE BldNdOuts_SetOutParam(BldNd_OutList, p, p_AD, ErrStat, ErrMsg ) endif ! The following are valid only for BEMT/DBEMT - if (p_AD%WakeMod /= WakeMod_BEMT) then + if (p_AD%Wake_Mod /= WakeMod_BEMT) then InvalidOutput( BldNd_BEM_F_qs ) = .true. InvalidOutput( BldNd_BEM_k_qs ) = .true. InvalidOutput( BldNd_BEM_kp_qs ) = .true. @@ -1553,7 +1553,7 @@ SUBROUTINE BldNdOuts_SetOutParam(BldNd_OutList, p, p_AD, ErrStat, ErrMsg ) ! it's going to be very difficult to get the FVW states without rewriting a bunch of code - if (.not. p_AD%UA_Flag .or. p_AD%WakeMod == WakeMod_FVW) then ! also invalid if AFAeroMod is not 4,5,6 + if (.not. p_AD%UA_Flag .or. p_AD%Wake_Mod == WakeMod_FVW) then ! also invalid if AFAeroMod is not 4,5,6 InvalidOutput( BldNd_UA_x1 ) = .true. InvalidOutput( BldNd_UA_x2 ) = .true. InvalidOutput( BldNd_UA_x3 ) = .true. diff --git a/modules/aerodyn/src/AeroDyn_Driver.f90 b/modules/aerodyn/src/AeroDyn_Driver.f90 index fdc343249c..076f0a4bcc 100644 --- a/modules/aerodyn/src/AeroDyn_Driver.f90 +++ b/modules/aerodyn/src/AeroDyn_Driver.f90 @@ -52,24 +52,24 @@ program AeroDyn_Driver ! Init of time estimator t_global=0.0_DbKi t_final=dat%dvr%numSteps*dat%dvr%dt - if (dat%dvr%analysisType/=idAnalysisCombi) then + !if (dat%dvr%analysisType/=idAnalysisCombi) then call SimStatus_FirstTime( TiLstPrn, PrevClockTime, SimStrtTime, UsrTime2, t_global, t_final ) - endif + !endif ! One time loop do nt = 1, dat%dvr%numSteps call Dvr_TimeStep(nt, dat%dvr, dat%ADI, dat%FED, dat%errStat, dat%errMsg); call CheckError() ! Time update to screen t_global=nt*dat%dvr%dt - if (dat%dvr%analysisType/=idAnalysisCombi) then + !if (dat%dvr%analysisType/=idAnalysisCombi) then if (mod( nt + 1, 10 )==0) call SimStatus(TiLstPrn, PrevClockTime, t_global, t_final) - endif + !endif end do !nt=1,numSteps - if (dat%dvr%analysisType/=idAnalysisCombi) then + !if (dat%dvr%analysisType/=idAnalysisCombi) then ! display runtime to screen call RunTimes(StrtTime, UsrTime1, SimStrtTime, UsrTime2, t_global) - endif + !endif call Dvr_EndCase(dat%dvr, dat%ADI, dat%initialized, dat%errStat, dat%errMsg); call CheckError() diff --git a/modules/aerodyn/src/AeroDyn_Driver_Subs.f90 b/modules/aerodyn/src/AeroDyn_Driver_Subs.f90 index 39f451791b..405cdb4ddb 100644 --- a/modules/aerodyn/src/AeroDyn_Driver_Subs.f90 +++ b/modules/aerodyn/src/AeroDyn_Driver_Subs.f90 @@ -445,7 +445,7 @@ subroutine Init_ADI_ForDriver(iCase, ADI, dvr, FED, dt, errStat, errMsg) !call AD_Dvr_DestroyAeroDyn_Data (AD , errStat2, errMsg2); call SetErrStat(errStat2, errMsg2, errStat, errMsg, RoutineName) needInit=.true. endif - if (ADI%p%AD%WakeMod == WakeMod_FVW) then + if (ADI%p%AD%Wake_Mod == WakeMod_FVW) then call WrScr('[INFO] OLAF is used, AeroDyn will be re-initialized') needInit=.true. endif @@ -496,15 +496,15 @@ subroutine Init_ADI_ForDriver(iCase, ADI, dvr, FED, dt, errStat, errMsg) if (wt%projMod==-1)then !call WrScr('>>> Using HAWTprojection to determine projMod') if (wt%HAWTprojection) then - InitInp%AD%rotors(iWT)%AeroProjMod = APM_BEM_NoSweepPitchTwist ! default, with WithoutSweepPitchTwist + !InitInp%AD%rotors(iWT)%AeroProjMod = APM_BEM_NoSweepPitchTwist ! default, with WithoutSweepPitchTwist + InitInp%AD%rotors(iWT)%AeroProjMod = -1 ! We let the code decide based on BEM_Mod else InitInp%AD%rotors(iWT)%AeroProjMod = APM_LiftingLine endif else InitInp%AD%rotors(iWT)%AeroProjMod = wt%projMod endif - InitInp%AD%rotors(iWT)%AeroBEM_Mod = wt%BEM_Mod - !call WrScr(' Driver: projMod: '//trim(num2lstr(InitInp%AD%rotors(iWT)%AeroProjMod))//', BEM_Mod:'//trim(num2lstr(InitInp%AD%rotors(iWT)%AeroBEM_Mod))) + call WrScr(' Driver: projMod: '//trim(num2lstr(InitInp%AD%rotors(iWT)%AeroProjMod))) InitInp%AD%rotors(iWT)%HubPosition = y_ED%HubPtMotion%Position(:,1) InitInp%AD%rotors(iWT)%HubOrientation = y_ED%HubPtMotion%RefOrientation(:,:,1) InitInp%AD%rotors(iWT)%NacellePosition = y_ED%NacelleMotion%Position(:,1) @@ -1009,9 +1009,6 @@ subroutine Dvr_ReadInputFile(fileName, dvr, errStat, errMsg ) call ParseVar(FileInfo_In, CurLine, 'ProjMod'//sWT , wt%projMod , errStat2, errMsg2, unEc); if (errStat2==ErrID_Fatal) then wt%projMod = -1 - wt%BEM_Mod = -1 - else - call ParseVar(FileInfo_In, CurLine, 'BEM_Mod'//sWT , wt%BEM_Mod , errStat2, errMsg2, unEc); if(Failed()) return endif call ParseVar(FileInfo_In, CurLine, 'BasicHAWTFormat'//sWT , wt%basicHAWTFormat , errStat2, errMsg2, unEc); if(Failed()) return diff --git a/modules/aerodyn/src/AeroDyn_IO.f90 b/modules/aerodyn/src/AeroDyn_IO.f90 index 0674d11787..36926fa627 100644 --- a/modules/aerodyn/src/AeroDyn_IO.f90 +++ b/modules/aerodyn/src/AeroDyn_IO.f90 @@ -123,7 +123,7 @@ SUBROUTINE Calc_WriteOutput( p, p_AD, u, x, m, m_AD, y, OtherState, xd, indx, iR tmpHubMB = 0.0_ReKi ! Compute max radius and rotor speed - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then rmax = 0.0_ReKi do k=1,p%NumBlades do j=1,p%NumBlNds @@ -141,7 +141,7 @@ SUBROUTINE Calc_WriteOutput( p, p_AD, u, x, m, m_AD, y, OtherState, xd, indx, iR call Calc_WriteOutput_AD() ! need to call this before calling the BEMT vs FVW versions of outputs so that the integrated output quantities are known - if (p_AD%WakeMod /= WakeMod_FVW) then + if (p_AD%Wake_Mod /= WakeMod_FVW) then call Calc_WriteOutput_BEMT() else call Calc_WriteOutput_FVW() @@ -658,13 +658,28 @@ SUBROUTINE ParsePrimaryFileInfo( PriPath, InitInp, InputFile, RootName, NumBlade character(ErrMsgLen) :: ErrMsg_NoAllBldNdOuts integer(IntKi) :: CurLine !< current entry in FileInfo_In%Lines array real(ReKi) :: TmpRe5(5) !< temporary 8 number array for reading values in - + character(1024) :: sDummy !< temporary string + character(1024) :: tmpOutStr !< temporary string for writing to screen + logical :: wakeModProvided, frozenWakeProvided, skewModProvided, AFAeroModProvided, UAModProvided, isLegalComment, firstWarn !< Temporary for legacy purposes + logical :: AoA34_Missing + integer :: UAMod_Old + integer :: WakeMod_Old + integer :: AFAeroMod_Old + integer :: SkewMod_Old + logical :: FrozenWake_Old character(*), parameter :: RoutineName = 'ParsePrimaryFileInfo' + UAMod_Old = -1 + WakeMod_Old = -1 + AFAeroMod_Old = -1 + SkewMod_Old = -1 + FrozenWake_Old = .False. + ! Initialization ErrStat = ErrId_None ErrMsg = "" UnEc = -1 ! Echo file unit. >0 when used + firstWarn=.False. CALL AllocAry( InputFileData%OutList, MaxOutPts, "Outlist", ErrStat2, ErrMsg2 ) @@ -702,12 +717,24 @@ SUBROUTINE ParsePrimaryFileInfo( PriPath, InitInp, InputFile, RootName, NumBlade ! DTAero - Time interval for aerodynamic calculations {or default} (s): call ParseVarWDefault ( FileInfo_In, CurLine, "DTAero", InputFileData%DTAero, interval, ErrStat2, ErrMsg2, UnEc ) if (Failed()) return - ! WakeMod - Type of wake/induction model (switch) {0=none, 1=BEMT, 2=DBEMT, 3=OLAF} [WakeMod cannot be 2 or 3 when linearizing] - call ParseVar( FileInfo_In, CurLine, "WakeMod", InputFileData%WakeMod, ErrStat2, ErrMsg2, UnEc ) - if (Failed()) return - ! AFAeroMod - Type of blade airfoil aerodynamics model (switch) {1=steady model, 2=Beddoes-Leishman unsteady model} [AFAeroMod must be 1 when linearizing] - call ParseVar( FileInfo_In, CurLine, "AFAeroMod", InputFileData%AFAeroMod, ErrStat2, ErrMsg2, UnEc ) - if (Failed()) return + ! WakeMod - LEGACY + call ParseVar( FileInfo_In, CurLine, "WakeMod", WakeMod_Old, ErrStat2, ErrMsg2, UnEc) + wakeModProvided = legacyInputPresent('WakeMod', CurLine, ErrStat2, ErrMsg2, 'Wake_Mod=0 (WakeMod=0), Wake_Mod=1 (WakeMod=1), DBEMT_Mod>0 (WakeMod=2), Wake_Mod=3 (WakeMod=3)') + ! Wake_Mod- Type of wake/induction model (switch) {0=none, 1=BEMT, 2=TBD, 3=OLAF} + call ParseVar( FileInfo_In, CurLine, "Wake_Mod", InputFileData%Wake_Mod, ErrStat2, ErrMsg2, UnEc ) + if (newInputMissing('Wake_Mod', CurLine, errStat2, errMsg2)) then + call WrScr(' Setting Wake_Mod to 1 (BEM active) as the input is Missing (typical behavior).') + InputFileData%Wake_Mod = WakeMod_BEMT + else + if (wakeModProvided) then + call LegacyAbort('Cannot have both Wake_Mod and WakeMod in the input file'); return + endif + endif + + + ! AFAeroMod - Type of blade airfoil aerodynamics model (switch) {1=steady model, 2=Beddoes-Leishman unsteady model} [AFAeroMod must be 1 when linearizing] + call ParseVar( FileInfo_In, CurLine, "AFAeroMod", AFAeroMod_Old, ErrStat2, ErrMsg2, UnEc ) + AFAeroModProvided = legacyInputPresent('AFAeroMod', CurLine, ErrStat2, ErrMsg2, 'UA_Mod=0 (AFAeroMod=1) or UA_Mod>1 (AFAeroMod=2)') ! TwrPotent - Type of tower influence on wind based on potential flow around the tower (switch) {0=none, 1=baseline potential flow, 2=potential flow with Bak correction} call ParseVar( FileInfo_In, CurLine, "TwrPotent", InputFileData%TwrPotent, ErrStat2, ErrMsg2, UnEc ) if (Failed()) return @@ -717,9 +744,9 @@ SUBROUTINE ParsePrimaryFileInfo( PriPath, InitInp, InputFile, RootName, NumBlade ! TwrAero - Calculate tower aerodynamic loads? (flag) call ParseVar( FileInfo_In, CurLine, "TwrAero", InputFileData%TwrAero, ErrStat2, ErrMsg2, UnEc ) if (Failed()) return - ! FrozenWake - Assume frozen wake during linearization? (flag) [used only when WakeMod=1 and when linearizing] - call ParseVar( FileInfo_In, CurLine, "FrozenWake", InputFileData%FrozenWake, ErrStat2, ErrMsg2, UnEc ) - if (Failed()) return + ! FrozenWake - Assume frozen wake during linearization? (flag) [used only when WakeMod=1 and when linearizing] + call ParseVar( FileInfo_In, CurLine, "FrozenWake", FrozenWake_Old, ErrStat2, ErrMsg2, UnEc ) + frozenWakeProvided = legacyInputPresent('FrozenWake', Curline, ErrStat2, ErrMsg2, 'DBEMTMod=-1 (FrozenWake=True) or DBEMTMod>-1 (FrozenWake=False)') ! CavitCheck - Perform cavitation check? (flag) [AFAeroMod must be 1 when CavitCheck=true] call ParseVar( FileInfo_In, CurLine, "CavitCheck", InputFileData%CavitCheck, ErrStat2, ErrMsg2, UnEc ) if (Failed()) return @@ -754,14 +781,55 @@ SUBROUTINE ParsePrimaryFileInfo( PriPath, InitInp, InputFile, RootName, NumBlade if (Failed()) return !====== Blade-Element/Momentum Theory Options ====================================================== [unused when WakeMod=0 or 3] - if ( InputFileData%Echo ) WRITE(UnEc, '(A)') FileInfo_In%Lines(CurLine) ! Write section break to echo - CurLine = CurLine + 1 - ! SkewMod - Type of skewed-wake correction model (switch) {1=uncoupled, 2=Pitt/Peters, 3=coupled} [unused when WakeMod=0 or 3] - call ParseVar( FileInfo_In, CurLine, "SkewMod", InputFileData%SkewMod, ErrStat2, ErrMsg2, UnEc ) - if (Failed()) return + call ParseCom (FileInfo_in, CurLine, sDummy, errStat2, errMsg2, UnEc, isLegalComment); if (Failed()) return + + ! BEM_Mod + call ParseVar( FileInfo_In, CurLine, "BEM_Mod", InputFileData%BEM_Mod, ErrStat2, ErrMsg2, UnEc ) + if (newInputMissing('BEM_Mod', CurLine, errStat2, errMsg2)) then + call WrScr(' Setting BEM_Mod to 1 (NoPitchSweepPitch) as the input is Missing (legacy behavior).') + InputFileData%BEM_Mod = BEMMod_2D + else + call ParseCom (FileInfo_in, CurLine, sDummy, errStat2, errMsg2, UnEc, isLegalComment); if (Failed()) return + endif + + ! SkewMod Legacy + call ParseVar( FileInfo_In, CurLine, "SkewMod", SkewMod_Old, ErrStat2, ErrMsg2, UnEc ) + skewModProvided = legacyInputPresent('SkewMod', CurLine, ErrStat2, ErrMsg2, 'Skew_Mod=-1 (SkewMod=0), Skew_Mod=0 (SkewMod=1), Skew_Mod=1 (SkewMod>=2)') + ! Skew_Mod- Select skew model {0: No skew model at all, -1:Throw away non-normal component for linearization, 1: Glauert skew model, } + call ParseVar( FileInfo_In, CurLine, "Skew_Mod", InputFileData%Skew_Mod, ErrStat2, ErrMsg2, UnEc ) + if (newInputMissing('Skew_Mod', CurLine, errStat2, errMsg2)) then + call WrScr(' Setting Skew_Mod to 1 (skew active) as the input is Missing (typical behavior).') + InputFileData%Skew_Mod = Skew_Mod_Active + else + if (skewModProvided) then + call LegacyAbort('Cannot have both Skew_Mod and SkewMod in the input file'); return + endif + endif + + + ! SkewMomCorr - Turn the skew momentum correction on or off [used only when SkewMod=1] + call ParseVar( FileInfo_In, CurLine, "SkewMomCorr", InputFileData%SkewMomCorr, ErrStat2, ErrMsg2, UnEc ) + if (newInputMissing('SkewMomCorr', CurLine, errStat2, errMsg2)) then + call WrScr(' Setting SkewMomCorr to False as the input is Missing (legacy behavior).') + InputFileData%SkewMomCorr = .False. + endif + + ! SkewRedistr_Mod - Type of skewed-wake correction model (switch) {0: no redistribution, 1=Glauert/Pitt/Peters, 2=Vortex Cylinder} [unsed only when SkewMod=1] + call ParseVarWDefault( FileInfo_In, CurLine, "SkewRedistr_Mod", InputFileData%SkewRedistr_Mod, 1, ErrStat2, ErrMsg2, UnEc ) + if (newInputMissing('SkewRedistr_Mod', CurLine, errStat2, errMsg2)) then + call WrScr(' Setting SkewRedistr_Mod to 1 as the input is Missing (legacy behavior).') + InputFileData%SkewRedistr_Mod = 1 + endif + ! SkewModFactor - Constant used in Pitt/Peters skewed wake model {or "default" is 15/32*pi} (-) [used only when SkewMod=2; unused when WakeMod=0 or 3] call ParseVarWDefault( FileInfo_In, CurLine, "SkewModFactor", InputFileData%SkewModFactor, (15.0_ReKi * pi / 32.0_ReKi), ErrStat2, ErrMsg2, UnEc ) - if (Failed()) return + if( legacyInputPresent('SkewModFactor', CurLine, ErrStat2, ErrMsg2, 'Rename this parameter to SkewRedistrFactor')) then + ! pass + else + call ParseVarWDefault( FileInfo_In, CurLine, "SkewRedistrFactor", InputFileData%SkewModFactor, (15.0_ReKi * pi / 32.0_ReKi), ErrStat2, ErrMsg2, UnEc ); if (Failed()) return + call ParseCom (FileInfo_in, CurLine, sDummy, errStat2, errMsg2, UnEc, isLegalComment); if (Failed()) return + endif + ! TipLoss - Use the Prandtl tip-loss model? (flag) [unused when WakeMod=0 or 3] call ParseVar( FileInfo_In, CurLine, "TipLoss", InputFileData%TipLoss, ErrStat2, ErrMsg2, UnEc ) if (Failed()) return @@ -787,11 +855,22 @@ SUBROUTINE ParsePrimaryFileInfo( PriPath, InitInp, InputFile, RootName, NumBlade ! MaxIter - Maximum number of iteration steps (-) [unused when WakeMod=0] call ParseVar( FileInfo_In, CurLine, "MaxIter", InputFileData%MaxIter, ErrStat2, ErrMsg2, UnEc ) if (Failed()) return + ! --- Shear + call ParseCom (FileInfo_in, CurLine, sDummy, errStat2, errMsg2, UnEc, isLegalComment); if (Failed()) return + call ParseVar( FileInfo_In, CurLine, "SectAvg" , InputFileData%SectAvg, ErrStat2, ErrMsg2, UnEc ); + if (newInputMissing('SectAvg', CurLine, errStat2, errMsg2)) then + call WrScr(' Setting SectAvg to False as the input is Missing (legacy behavior).') + InputFileData%SectAvg = .false. + else + call ParseVarWDefault( FileInfo_In, CurLine, "SectAvgWeighting", InputFileData%SA_Weighting, 1 , ErrStat2, ErrMsg2, UnEc ); if (Failed()) return + call ParseVarWDefault( FileInfo_In, CurLine, "SectAvgNPoints" , InputFileData%SA_nPerSec, 5 , ErrStat2, ErrMsg2, UnEc ); if (Failed()) return + call ParseVarWDefault( FileInfo_In, CurLine, "SectAvgPsiBwd" , InputFileData%SA_PsiBwd, -60._ReKi, ErrStat2, ErrMsg2, UnEc ); if (Failed()) return + call ParseVarWDefault( FileInfo_In, CurLine, "SectAvgPsiFwd" , InputFileData%SA_PsiFwd, 60._ReKi, ErrStat2, ErrMsg2, UnEc ); if (Failed()) return + call ParseCom (FileInfo_in, CurLine, sDummy, errStat2, errMsg2, UnEc, isLegalComment); if (Failed()) return + endif !====== Dynamic Blade-Element/Momentum Theory Options ============================================== [used only when WakeMod=2] - if ( InputFileData%Echo ) WRITE(UnEc, '(A)') FileInfo_In%Lines(CurLine) ! Write section break to echo - CurLine = CurLine + 1 - ! DBEMT_Mod - Type of dynamic BEMT (DBEMT) model {1=constant tau1, 2=time-dependent tau1} (-) [used only when WakeMod=2] + ! DBEMT_Mod - Type of dynamic BEMT (DBEMT) model {1=constant tau1, 2=time-dependent tau1} (-) [used only when WakeMod=2] call ParseVar( FileInfo_In, CurLine, "DBEMT_Mod", InputFileData%DBEMT_Mod, ErrStat2, ErrMsg2, UnEc ) if (Failed()) return ! tau1_const - Time constant for DBEMT (s) [used only when WakeMod=2 and DBEMT_Mod=1] @@ -807,11 +886,28 @@ SUBROUTINE ParsePrimaryFileInfo( PriPath, InitInp, InputFile, RootName, NumBlade IF ( PathIsRelative( InputFileData%FVWFileName ) ) InputFileData%FVWFileName = TRIM(PriPath)//TRIM(InputFileData%FVWFileName) !====== Beddoes-Leishman Unsteady Airfoil Aerodynamics Options ===================================== [used only when AFAeroMod=2] - if ( InputFileData%Echo ) WRITE(UnEc, '(A)') FileInfo_In%Lines(CurLine) ! Write section break to echo - CurLine = CurLine + 1 - ! UAMod - Unsteady Aero Model Switch (switch) {1=Baseline model (Original), 2=Gonzalez's variant (changes in Cn,Cc,Cm), 3=Minnema/Pierce variant (changes in Cc and Cm)} [used only when AFAeroMod=2] - call ParseVar( FileInfo_In, CurLine, "UAMod", InputFileData%UAMod, ErrStat2, ErrMsg2, UnEc ) - if (Failed()) return + call ParseCom (FileInfo_in, CurLine, sDummy, errStat2, errMsg2, UnEc, isLegalComment); if (Failed()) return + ! AoA34 Sample the angle of attack (AoA) at the 3/4 chord or the AC point {default=True} [always used] + call ParseVar( FileInfo_In, CurLine, "AoA34", InputFileData%AoA34, ErrStat2, ErrMsg2, UnEc ) + AoA34_Missing = newInputMissing('AoA34', CurLine, errStat2, errMsg2) + ! UAMod (Legacy) + call ParseVar( FileInfo_In, CurLine, "UAMod", UAMod_Old, ErrStat2, ErrMsg2, UnEc ) + UAModProvided = legacyInputPresent('UAMod', CurLine, ErrStat2, ErrMsg2, 'UA_Mod=0 (AFAeroMod=1), UA_Mod>1 (AFAeroMod=2 and UA_Mod=UAMod') + ! UA_Mod - Unsteady Aero Model Switch (switch) {0=Quasi-steady (no UA), 2=Gonzalez's variant (changes in Cn,Cc,Cm), 3=Minnema/Pierce variant (changes in Cc and Cm)} + call ParseVar( FileInfo_In, CurLine, "UA_Mod", InputFileData%UA_Mod, ErrStat2, ErrMsg2, UnEc ) + if (newInputMissing('UA_Mod', CurLine, errStat2, errMsg2)) then + ! We'll deal with it when we deal with AFAeroMod + InputFileData%UA_Mod = UAMod_Old + if (.not. UAModProvided) then + call LegacyAbort('Need to provide either UA_Mod or UAMod in the input file'); return + endif + else + if (UAModProvided) then + call LegacyAbort('Cannot have both UA_Mod and UAMod in the input file'); return + endif + endif + + ! FLookup - Flag to indicate whether a lookup for f' will be calculated (TRUE) or whether best-fit exponential equations will be used (FALSE); if FALSE S1-S4 must be provided in airfoil input files (flag) [used only when AFAeroMod=2] call ParseVar( FileInfo_In, CurLine, "FLookup", InputFileData%FLookup, ErrStat2, ErrMsg2, UnEc ) if (Failed()) return @@ -987,6 +1083,62 @@ SUBROUTINE ParsePrimaryFileInfo( PriPath, InitInp, InputFile, RootName, NumBlade call ReadOutputListFromFileInfo( FileInfo_In, CurLine, InputFileData%OutList, InputFileData%NumOuts, ErrStat2, ErrMsg2, UnEc ) if (Failed()) return; + !====== Legacy logic to match old and new input files ================================================ + ! NOTE: remove me in future release + if (frozenWakeProvided) then + if (FrozenWake_Old) then + call WrScr('> FrozenWake=True -> Setting DBEMT_Mod=-1') + InputFileData%DBEMT_Mod = DBEMT_frozen + else + call WrScr('> FrozenWake=False -> Not changing DBEMT_Mod') + endif + endif + if (wakeModProvided) then + InputFileData%Wake_Mod = WakeMod_Old + if (WakeMod_Old==1) then + call WrScr('> WakeMod=1 -> Setting DBEMT_Mod=0') + ! Turn off DBEMT + InputFileData%DBEMT_Mod=DBEMT_none + else if (WakeMod_Old==2) then + call WrScr('> WakeMod=2 -> Setting Wake_Mod=1 (BEMT) (DBEMT_Mod needs to be >0)') + InputFileData%Wake_Mod = WakeMod_BEMT + if (InputFileData%DBEMT_Mod < DBEMT_none) then + call LegacyAbort('DBEMT should be >0 when using legacy input WakeMod=2'); return + endif + endif + endif + if (AFAeroModProvided) then + if (AFAeroMod_Old==1) then + call WrScr('> AFAeroMod=1 -> Setting UA_Mod=0') + InputFileData%UA_Mod = UA_None + if (AoA34_Missing) then + call WrScr('> Setting AoA34 to False as the input is Missing and UA is turned off (legacy behavior).') + InputFileData%AoA34=.false. + endif + else if (AFAeroMod_Old==2) then + call WrScr('> AFAeroMod=2 -> Not changing DBEMT_Mod') + if (InputFileData%UA_Mod==0) then + call LegacyAbort('Cannot set UA_Mod=0 with legacy option AFAeroMod=2 (inconsistent behavior).'); return + else if (AoA34_Missing) then + call WrScr('> Setting AoA34 to True as the input is Missing and UA is turned on (legacy behavior).') + InputFileData%AoA34=.true. + endif + else + call LegacyAbort('AFAeroMod should be 1 or 2'); return + endif + endif + if (skewModProvided) then + if (SkewMod_Old==0) then + InputFileData%Skew_Mod = Skew_Mod_Orthogonal + else if (SkewMod_Old==1) then + InputFileData%Skew_Mod = Skew_Mod_None + else if (SkewMod_Old==2) then + InputFileData%Skew_Mod = Skew_Mod_Active + else + call LegacyAbort('Legacy option SkewMod is not 0, 1,2 which is not supported.'); return + endif + endif + !====== Nodal Outputs ============================================================================== ! In case there is something ill-formed in the additional nodal outputs section, we will simply ignore it. ! Expecting at least 5 more lines in the input file for this section @@ -1019,6 +1171,41 @@ SUBROUTINE ParsePrimaryFileInfo( PriPath, InitInp, InputFile, RootName, NumBlade ! Prevent segfault when no blades specified. All logic tests on BldNd_NumOuts at present. if (InputFileData%BldNd_BladesOut <= 0) InputFileData%BldNd_NumOuts = 0 + + !====== Print new and old inputs ===================================================================== + if (wakeModProvided .or. frozenWakeProvided .or. skewModProvided .or. AFAeroModProvided .or. UAModProvided) then + call printNewOldInputs() + endif + + !====== Advanced Options ============================================================================= + if ((CurLine) >= size(FileInfo_In%Lines)) RETURN + + call WrScr(' - Reading advanced options for AeroDyn') + do CurLine= CurLine, size(FileInfo_In%Lines) + sDummy = FileInfo_In%Lines(CurLine) + call Conv2UC(sDummy) ! to uppercase + if (index(sDummy, '!') == 1 .or. index(sDummy, '=') == 1 .or. index(sDummy, '#') == 1 .or. index(sDummy, '---') == 1) then + ! pass comment lines + elseif (index(sDummy, 'SECTAVG')>1) then + read(sDummy, '(L1)') InputFileData%SectAvg + write(tmpOutStr,*) ' >>> SectAvg ',InputFileData%SectAvg + elseif (index(sDummy, 'SA_PSIBWD')>1) then + read(sDummy, *) InputFileData%SA_PsiBwd + write(tmpOutStr,*) ' >>> SA_PsiBwd ',InputFileData%SA_PsiBwd + elseif (index(sDummy, 'SA_PSIFWD')>1) then + read(sDummy, *) InputFileData%SA_PsiFwd + write(tmpOutStr,*) ' >>> SA_PsiFwd ',InputFileData%SA_PsiFwd + elseif (index(sDummy, 'SA_NPERSEC')>1) then + read(sDummy, *) InputFileData%SA_nPerSec + write(tmpOutStr,*) ' >>> SA_nPerSec ',InputFileData%SA_nPerSec + else + write(tmpOutStr,*) '[WARN] AeroDyn Line ignored: '//trim(sDummy) + endif + call WrScr(trim(tmpOutStr)) + enddo + + + RETURN CONTAINS !------------------------------------------------------------------------------------------------- @@ -1030,22 +1217,107 @@ logical function Failed() end function Failed logical function FailedNodal() ErrMsg_NoAllBldNdOuts='AD15 Nodal Outputs: Nodal output section of AeroDyn input file not found or improperly formatted. Skipping nodal outputs.' - FailedNodal = ErrStat2 >= AbortErrLev + ! TODO Use and ErrID_Fatal here + FailedNodal = ErrStat2 >= AbortErrLev if ( FailedNodal ) then InputFileData%BldNd_BladesOut = 0 InputFileData%BldNd_NumOuts = 0 call wrscr( trim(ErrMsg_NoAllBldNdOuts) ) + call printNewOldInputs() endif end function FailedNodal subroutine LegacyWarning(Message) character(len=*), intent(in) :: Message - call WrScr('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') - call WrScr('Warning: the AeroDyn input file is not at the latest format!' ) - call WrScr(' Visit: https://openfast.readthedocs.io/en/dev/source/user/api_change.html') + if (.not.FirstWarn) then + call WrScr('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') + call WrScr('[WARN] The AeroDyn input file is not at the latest format!' ) + call WrScr(' Visit: https://openfast.readthedocs.io/en/dev/source/user/api_change.html') + call WrScr('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') + FirstWarn=.true. + endif call WrScr('> Issue: '//trim(Message)) - call WrScr('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') end subroutine LegacyWarning !------------------------------------------------------------------------------------------------- + subroutine LegacyAbort(Message) + character(len=*), intent(in) :: Message + call SetErrStat( ErrID_Fatal, Message, ErrStat, ErrMsg, 'ParsePrimaryFileInfo' ) + end subroutine LegacyAbort + !------------------------------------------------------------------------------------------------- + logical function legacyInputPresent(varName, iLine, errStat, errMsg, varNameSubs) + character(len=*), intent(in ) :: varName !< Variable being read + integer(IntKi), intent(in ) :: iLine !< Line number + integer(IntKi), intent(inout) :: errStat !< Error status + character(ErrMsgLen), intent(inout) :: errMsg !< Error message + character(len=*), optional, intent(in ) :: varNameSubs !< Substituted variable + legacyInputPresent = errStat == ErrID_None + if (legacyInputPresent) then + if (present(varNameSubs)) then + call LegacyWarning(trim(varName)//' has now been removed.'//NewLine//' Use: '//trim(varNameSubs)//'.') + else + call LegacyWarning(trim(varName)//' has now been removed.') + endif + else + ! We are actually happy, this input should indeed not be present. + endif + ! We erase the error no matter what + errStat = ErrID_None + errMsg = '' + end function legacyInputPresent + !------------------------------------------------------------------------------------------------- + logical function newInputMissing(varName, iLine, errStat, errMsg, varNameSubs) + character(len=*), intent(in ) :: varName !< Variable being read + integer(IntKi), intent(in ) :: iLine !< Line number + integer(IntKi), intent(inout) :: errStat !< Error status + character(ErrMsgLen), intent(inout) :: errMsg !< Error message + character(len=*), optional, intent(in ) :: varNameSubs !< Substituted variable + newInputMissing = errStat == ErrID_Fatal + if (newInputMissing) then + call LegacyWarning(trim(varName)//' should be present on line '//trim(num2lstr(iLine))//'.') + else + ! We are happy + endif + ! We erase the error + errStat = ErrID_None + errMsg = '' + end function newInputMissing + + !------------------------------------------------------------------------------------------------- + subroutine printNewOldInputs() + character(1024) :: tmpStr + ! Temporary HACK, for WakeMod=10, 11 or 12 use AeroProjMod 2 (will trigger PolarBEM) + if (InputFileData%Wake_Mod==10) then + call WrScr('[WARN] Wake_Mod=10 is a temporary hack. Setting BEM_Mod to 0') + InputFileData%BEM_Mod = 0 + elseif (InputFileData%Wake_Mod==11) then + call WrScr('[WARN] Wake_Mod=11 is a temporary hack. Setting BEM_Mod to 2') + InputFileData%BEM_Mod = 2 + elseif (InputFileData%Wake_Mod==12) then + call WrScr('[WARN] Wake_Mod=12 is a temporary hack. Setting BEM_Mod to 2') + InputFileData%BEM_Mod = 2 + endif + !====== Summary of new AeroDyn options =============================================================== + ! NOTE: remove me in future release + call WrScr('-------------- New AeroDyn inputs (with new meaning):') + write (tmpStr,'(A20,I0)') 'Wake_Mod: ' , InputFileData%Wake_Mod; call WrScr(trim(tmpStr)) + write (tmpStr,'(A20,I0)') 'BEM_Mod: ' , InputFileData%BEM_Mod; call WrScr(trim(tmpStr)) + write (tmpStr,'(A20,L1)') 'SectAvg: ' , InputFileData%SectAvg; call WrScr(trim(tmpStr)) + write (tmpStr,'(A20,I0)') 'SectAvgWeighting: ', InputFileData%SA_Weighting; call WrScr(trim(tmpStr)) + write (tmpStr,'(A20,I0)') 'SectAvgNPoints: ', InputFileData%SA_nPerSec; call WrScr(trim(tmpStr)) + write (tmpStr,'(A20,I0)') 'DBEMT_Mod:' , InputFileData%DBEMT_Mod; call WrScr(trim(tmpStr)) + write (tmpStr,'(A20,I0)') 'Skew_Mod: ' , InputFileData%Skew_Mod; call WrScr(trim(tmpStr)) + write (tmpStr,'(A20,L1)') 'SkewMomCorr:' , InputFileData%SkewMomCorr; call WrScr(trim(tmpStr)) + write (tmpStr,'(A20,I0)') 'SkewRedistr_Mod:' , InputFileData%SkewRedistr_Mod; call WrScr(trim(tmpStr)) + write (tmpStr,'(A20,L1)') 'AoA34: ' , InputFileData%AoA34; call WrScr(trim(tmpStr)) + write (tmpStr,'(A20,I0)') 'UA_Mod: ' , InputFileData%UA_Mod; call WrScr(trim(tmpStr)) + call WrScr('-------------- Old AeroDyn inputs:') + write (tmpStr,'(A20,I0)') 'WakeMod: ', WakeMod_Old; call WrScr(trim(tmpStr)) + write (tmpStr,'(A20,I0)') 'SkewMod: ', SkewMod_Old; call WrScr(trim(tmpStr)) + write (tmpStr,'(A20,I0)') 'AFAeroMod:', AFAeroMod_Old; call WrScr(trim(tmpStr)) + write (tmpStr,'(A20,L1)') 'FrozenWake:', FrozenWake_Old; call WrScr(trim(tmpStr)) + write (tmpStr,'(A20,I0)') 'UAMod: ', UAMod_Old; call WrScr(trim(tmpStr)) + call WrScr('------------------------------------------------------') + end subroutine printNewOldInputs + END SUBROUTINE ParsePrimaryFileInfo !---------------------------------------------------------------------------------------------------------------------------------- SUBROUTINE ReadBladeInputs ( ADBlFile, BladeKInputFileData, AeroProjMod, UnEc, ErrStat, ErrMsg ) @@ -1339,11 +1611,9 @@ SUBROUTINE AD_PrintSum( InputFileData, p, p_AD, u, y, ErrStat, ErrMsg ) WRITE (UnSu,'(/,A)') '====== General Options ============================================================================' ! WakeMod - select case (p_AD%WakeMod) + select case (p_AD%Wake_Mod) case (WakeMod_BEMT) Msg = 'Blade-Element/Momentum Theory' - case (WakeMod_DBEMT) - Msg = 'Dynamic Blade-Element/Momentum Theory' case (WakeMod_FVW) Msg = 'Free Vortex Wake Theory' case (WakeMod_None) @@ -1351,20 +1621,7 @@ SUBROUTINE AD_PrintSum( InputFileData, p, p_AD, u, y, ErrStat, ErrMsg ) case default Msg = 'unknown' end select - WRITE (UnSu,Ec_IntFrmt) p_AD%WakeMod, 'WakeMod', 'Type of wake/induction model: '//TRIM(Msg) - - - ! AFAeroMod - select case (InputFileData%AFAeroMod) - case (AFAeroMod_BL_unsteady) - Msg = 'Beddoes-Leishman unsteady model' - case (AFAeroMod_steady) - Msg = 'steady' - case default - Msg = 'unknown' - end select - WRITE (UnSu,Ec_IntFrmt) InputFileData%AFAeroMod, 'AFAeroMod', 'Type of blade airfoil aerodynamics model: '//TRIM(Msg) - + WRITE (UnSu,Ec_IntFrmt) p_AD%Wake_Mod, 'WakeMod', 'Type of wake/induction model: '//TRIM(Msg) ! TwrPotent select case (p%TwrPotent) @@ -1419,21 +1676,21 @@ SUBROUTINE AD_PrintSum( InputFileData, p, p_AD, u, y, ErrStat, ErrMsg ) WRITE (UnSu,Ec_LgFrmt) p%Buoyancy, 'Buoyancy', 'Include buoyancy effects? '//TRIM(Msg) - if (p_AD%WakeMod/=WakeMod_none) then + if (p_AD%Wake_Mod/=WakeMod_none) then WRITE (UnSu,'(A)') '====== Blade-Element/Momentum Theory Options ======================================================' ! SkewMod - select case (InputFileData%SkewMod) - case (SkewMod_Orthogonal) + select case (InputFileData%Skew_Mod) + case (Skew_Mod_Orthogonal) Msg = 'orthogonal' - case (SkewMod_Uncoupled) - Msg = 'uncoupled' - case (SkewMod_PittPeters) - Msg = 'Pitt/Peters' + case (Skew_Mod_None) + Msg = 'no correction' + case (Skew_Mod_Active) + Msg = 'active' case default Msg = 'unknown' end select - WRITE (UnSu,Ec_IntFrmt) InputFileData%SkewMod, 'SkewMod', 'Type of skewed-wake correction model: '//TRIM(Msg) + WRITE (UnSu,Ec_IntFrmt) InputFileData%Skew_Mod, 'Skew_Mod', 'Skewed-wake correction model: '//TRIM(Msg) ! TipLoss @@ -1485,63 +1742,66 @@ SUBROUTINE AD_PrintSum( InputFileData, p, p_AD, u, y, ErrStat, ErrMsg ) ! MaxIter - if (p_AD%WakeMod == WakeMod_DBEMT) then - select case (InputFileData%DBEMT_Mod) - case (DBEMT_tauConst) - Msg = 'constant tau1' - case (DBEMT_tauVaries) - Msg = 'time-dependent tau1' - case (DBEMT_cont_tauConst) - Msg = 'continuous formulation with constant tau1' - case default - Msg = 'unknown' - end select - - WRITE (UnSu,Ec_IntFrmt) InputFileData%DBEMT_Mod, 'DBEMT_Mod', 'Type of dynamic BEMT (DBEMT) model: '//TRIM(Msg) + select case (InputFileData%DBEMT_Mod) + case (DBEMT_frozen) + Msg = 'frozen-wake' + case (DBEMT_none) + Msg = 'quasi-steady' + case (DBEMT_tauConst) + Msg = 'dynamic - constant tau1' + case (DBEMT_tauVaries) + Msg = 'dynamic - time-dependent tau1' + case (DBEMT_cont_tauConst) + Msg = 'dynamic - continuous formulation with constant tau1' + case default + Msg = 'unknown' + end select + + WRITE (UnSu,Ec_IntFrmt) InputFileData%DBEMT_Mod, 'DBEMT_Mod', 'Type of dynamic BEMT (DBEMT) model: '//TRIM(Msg) - if (InputFileData%DBEMT_Mod==DBEMT_tauConst) & - WRITE (UnSu,Ec_ReFrmt) InputFileData%tau1_const, 'tau1_const', 'Time constant for DBEMT (s)' + if (InputFileData%DBEMT_Mod==DBEMT_tauConst) & + WRITE (UnSu,Ec_ReFrmt) InputFileData%tau1_const, 'tau1_const', 'Time constant for DBEMT (s)' - end if end if - if (InputFileData%AFAeroMod==AFAeroMod_BL_unsteady) then - WRITE (UnSu,'(A)') '====== Beddoes-Leishman Unsteady Airfoil Aerodynamics Options =====================================' - - ! UAMod - select case (InputFileData%UAMod) - case (UA_Baseline) - Msg = 'baseline model (original)' - case (UA_Gonzalez) - Msg = "Gonzalez's variant (changes in Cn, Cc, and Cm)" - case (UA_MinnemaPierce) - Msg = 'Minnema/Pierce variant (changes in Cc and Cm)' - !case (4) - ! Msg = 'DYSTOOL' - case (UA_HGM) - Msg = 'HGM (continuous state)' - case (UA_HGMV) - Msg = 'HGMV (continuous state + vortex)' - case (UA_OYE) - Msg = 'Stieg Oye dynamic stall model' - case (UA_BV) - Msg = 'Boeing-Vertol dynamic stall model (e.g. used in CACTUS)' - case default - Msg = 'unknown' - end select - WRITE (UnSu,Ec_IntFrmt) InputFileData%UAMod, 'UAMod', 'Unsteady Aero Model: '//TRIM(Msg) + WRITE (UnSu,'(A)') '======================== Unsteady Airfoil Aerodynamics Options =====================================' - - ! FLookup - if (InputFileData%FLookup) then - Msg = 'Yes' - else - Msg = 'No, use best-fit exponential equations instead' - end if - WRITE (UnSu,Ec_LgFrmt) InputFileData%FLookup, 'FLookup', "Use a lookup for f'? "//TRIM(Msg) + ! UAMod + select case (InputFileData%UA_Mod) + case (UA_None) + Msg = 'none (quasi-steady airfoil aerodynamics)' + case (UA_Baseline) + Msg = 'baseline model (original)' + case (UA_Gonzalez) + Msg = "Gonzalez's variant (changes in Cn, Cc, and Cm)" + case (UA_MinnemaPierce) + Msg = 'Minnema/Pierce variant (changes in Cc and Cm)' + !case (4) + ! Msg = 'DYSTOOL' + case (UA_HGM) + Msg = 'HGM (continuous state)' + case (UA_HGMV) + Msg = 'HGMV (continuous state + vortex)' + case (UA_OYE) + Msg = 'Stieg Oye dynamic stall model' + case (UA_BV) + Msg = 'Boeing-Vertol dynamic stall model (e.g. used in CACTUS)' + case default + Msg = 'unknown' + end select + WRITE (UnSu,Ec_IntFrmt) InputFileData%UA_Mod, 'UA_Mod', 'Unsteady Aero Model: '//TRIM(Msg) + + + ! FLookup + if (InputFileData%FLookup) then + Msg = 'Yes' + else + Msg = 'No, use best-fit exponential equations instead' + end if + WRITE (UnSu,Ec_LgFrmt) InputFileData%FLookup, 'FLookup', "Use a lookup for f'? "//TRIM(Msg) + - end if WRITE (UnSu,'(A)') '====== Outputs ====================================================================================' @@ -1656,7 +1916,7 @@ SUBROUTINE SetOutParam(OutList, p, p_AD, ErrStat, ErrMsg ) end if - if (p_AD%WakeMod /= WakeMod_DBEMT) then + if (p%DBEMT_Mod > DBEMT_none) then InvalidOutput( DBEMTau1 ) = .true. end if @@ -2145,19 +2405,14 @@ subroutine calcCantAngle(f, xi,stencilSize,n,cantAngle) implicit none integer(IntKi), intent(in) :: stencilSize, n integer(IntKi) :: i, j - integer(IntKi) :: sortInd(n) integer(IntKi) :: info real(ReKi), intent(in) :: f(n), xi(n) real(ReKi) :: cx(stencilSize), cf(stencilSize), xiIn(stencilSize) - real(ReKi) :: fIn(stencilSize), cPrime(n), fPrime(n), xiAbs(n) + real(ReKi) :: fIn(stencilSize), cPrime(n), fPrime(n) real(ReKi), intent(inout) :: cantAngle(n) - real(ReKi), parameter :: ThisTol = 1e-6 do i = 1,size(xi) - xiAbs = abs(xi-xi(i)) - call hpsort_eps_epw (n, xiAbs, sortInd, ThisTol) - if (i.eq.1) then fIn = f(1:stencilSize) xiIn = xi(1:stencilSize) @@ -2171,7 +2426,7 @@ subroutine calcCantAngle(f, xi,stencilSize,n,cantAngle) call differ_stencil ( xi(i), 1, 2, xiIn, cx, info ) !call differ_stencil ( xi(i), 1, 2, fIn, cf, info ) if (info /= 0) then - print*,'Cant Calc failed at i=',i + call WrScr('Cant Calc failed at i='//trim(Num2LStr(i))) else cPrime(i) = 0.0 fPrime(i) = 0.0 @@ -2394,147 +2649,5 @@ subroutine r8vm_sl ( n, a, b, x, job, info ) return end subroutine r8vm_sl -! - ! Copyright (C) 2010-2016 Samuel Ponce', Roxana Margine, Carla Verdi, Feliciano Giustino - ! Copyright (C) 2007-2009 Jesse Noffsinger, Brad Malone, Feliciano Giustino - ! - ! This file is distributed under the terms of the GNU General Public - ! License. See the file `LICENSE' in the root directory of the - ! present distribution, or http://www.gnu.org/copyleft.gpl.txt . - ! - ! Adapted from flib/hpsort_eps - !--------------------------------------------------------------------- - subroutine hpsort_eps_epw (n, ra, ind, eps) - !--------------------------------------------------------------------- - ! sort an array ra(1:n) into ascending order using heapsort algorithm, - ! and considering two elements being equal if their values differ - ! for less than "eps". - ! n is input, ra is replaced on output by its sorted rearrangement. - ! create an index table (ind) by making an exchange in the index array - ! whenever an exchange is made on the sorted data array (ra). - ! in case of equal values in the data array (ra) the values in the - ! index array (ind) are used to order the entries. - ! if on input ind(1) = 0 then indices are initialized in the routine, - ! if on input ind(1) != 0 then indices are assumed to have been - ! initialized before entering the routine and these - ! indices are carried around during the sorting process - ! - ! no work space needed ! - ! free us from machine-dependent sorting-routines ! - ! - ! adapted from Numerical Recipes pg. 329 (new edition) - ! - !use kinds, ONLY : DP - implicit none - !-input/output variables - integer(IntKi), intent(in) :: n - real(ReKi), intent(in) :: eps - integer(IntKi) :: ind (n) - real(ReKi) :: ra (n) - !-local variables - integer(IntKi) :: i, ir, j, l, iind - real(ReKi) :: rra -! - ! initialize index array - IF (ind (1) .eq.0) then - DO i = 1, n - ind (i) = i - ENDDO - ENDIF - ! nothing to order - IF (n.lt.2) return - ! initialize indices for hiring and retirement-promotion phase - l = n / 2 + 1 - - ir = n - - sorting: do - - ! still in hiring phase - IF ( l .gt. 1 ) then - l = l - 1 - rra = ra (l) - iind = ind (l) - ! in retirement-promotion phase. - ELSE - ! clear a space at the end of the array - rra = ra (ir) - ! - iind = ind (ir) - ! retire the top of the heap into it - ra (ir) = ra (1) - ! - ind (ir) = ind (1) - ! decrease the size of the corporation - ir = ir - 1 - ! done with the last promotion - IF ( ir .eq. 1 ) then - ! the least competent worker at all ! - ra (1) = rra - ! - ind (1) = iind - exit sorting - ENDIF - ENDIF - ! wheter in hiring or promotion phase, we - i = l - ! set up to place rra in its proper level - j = l + l - ! - DO while ( j .le. ir ) - IF ( j .lt. ir ) then - ! compare to better underling - IF ( hslt( ra (j), ra (j + 1) ) ) then - j = j + 1 - !else if ( .not. hslt( ra (j+1), ra (j) ) ) then - ! this means ra(j) == ra(j+1) within tolerance - ! if (ind (j) .lt.ind (j + 1) ) j = j + 1 - ENDIF - ENDIF - ! demote rra - IF ( hslt( rra, ra (j) ) ) then - ra (i) = ra (j) - ind (i) = ind (j) - i = j - j = j + j - !else if ( .not. hslt ( ra(j) , rra ) ) then - !this means rra == ra(j) within tolerance - ! demote rra - ! if (iind.lt.ind (j) ) then - ! ra (i) = ra (j) - ! ind (i) = ind (j) - ! i = j - ! j = j + j - ! else - ! set j to terminate do-while loop - ! j = ir + 1 - ! endif - ! this is the right place for rra - ELSE - ! set j to terminate do-while loop - j = ir + 1 - ENDIF - ENDDO - ra (i) = rra - ind (i) = iind - - END DO sorting -contains - - ! internal function - ! compare two real number and return the result - - logical function hslt( a, b ) - REAL(ReKi) :: a, b - IF( abs(a-b) < eps ) then - hslt = .false. - ELSE - hslt = ( a < b ) - end if - end function hslt - - ! -end subroutine hpsort_eps_epw - !---------------------------------------------------------------------------------------------------------------------------------- END MODULE AeroDyn_IO diff --git a/modules/aerodyn/src/AeroDyn_Registry.txt b/modules/aerodyn/src/AeroDyn_Registry.txt index 6cfb51abd6..4161e04e2b 100644 --- a/modules/aerodyn/src/AeroDyn_Registry.txt +++ b/modules/aerodyn/src/AeroDyn_Registry.txt @@ -21,12 +21,9 @@ usefrom InflowWind.txt param AeroDyn/AD - IntKi ModelUnknown - -1 - "" - param ^ - IntKi WakeMod_none - 0 - "Wake model - none" - param ^ - IntKi WakeMod_BEMT - 1 - "Wake model - BEMT (blade elememnt momentum theory)" - -param ^ - IntKi WakeMod_DBEMT - 2 - "Wake model - DBEMT (dynamic elememnt momentum theory)" - +#param ^ - IntKi WakeMod_TODO - 2 - "Wake model - TBD" - param ^ - IntKi WakeMod_FVW - 3 - "Wake model - FVW (free vortex wake, OLAF)" - -param ^ - IntKi AFAeroMod_steady - 1 - "steady model" - -param ^ - IntKi AFAeroMod_BL_unsteady - 2 - "Beddoes-Leishman unsteady model" - - param ^ - IntKi TwrPotent_none - 0 - "no tower potential flow" - param ^ - IntKi TwrPotent_baseline - 1 - "baseline tower potential flow" - param ^ - IntKi TwrPotent_Bak - 2 - "tower potential flow with Bak correction" - @@ -35,6 +32,9 @@ param ^ - IntKi TwrShadow_none - 0 - "no tower s param ^ - IntKi TwrShadow_Powles - 1 - "Powles tower shadow model" - param ^ - IntKi TwrShadow_Eames - 2 - "Eames tower shadow model" - +param ^ - IntKi SA_Wgt_Uniform - 1 - "Sector average weighting - Uniform" - +#param ^ - IntKi SA_Wgt_Impulse - 1 - "Sector average weighting - Impulse" - + param ^ - IntKi TFinAero_none - 0 - "no tail fin aero" - param ^ - IntKi TFinAero_polar - 1 - "polar-based tail fin aerodynamics" - param ^ - IntKi TFinAero_USB - 2 - "unsteady slender body tail fin aerodynamics model" - @@ -95,7 +95,6 @@ typedef ^ RotInitInputType R8Ki BladeRootOrientation {:}{:}{:} - - "DCM referenc typedef ^ RotInitInputType R8Ki NacellePosition {3} - - "X-Y-Z reference position of nacelle" m typedef ^ RotInitInputType R8Ki NacelleOrientation {3}{3} - - "DCM reference orientation of nacelle" - typedef ^ RotInitInputType IntKi AeroProjMod - 1 - "Flag to switch between different projection models" - -typedef ^ RotInitInputType IntKi AeroBEM_Mod - -1 - "Flag to switch between different BEM Model" - typedef ^ RotInitInputType ReKi RotSpeed - - 0 "Rotor speed used when AeroDyn is computing aero maps" "rad/s" typedef ^ InitInputType RotInitInputType rotors {:} - - "Init Input Types for rotors" - @@ -176,12 +175,11 @@ typedef ^ RotInputFile TFinInputFileType TFin - - - "Input file typedef ^ AD_InputFile Logical Echo - - - "Echo input file to echo file" - typedef ^ AD_InputFile DbKi DTAero - - - "Time interval for aerodynamic calculations {or \"default\"}" s -typedef ^ AD_InputFile IntKi WakeMod - - - "Type of wake/induction model {0=none, 1=BEMT, 2=DBEMT, 3=FVW}" - -typedef ^ AD_InputFile IntKi AFAeroMod - - - "Type of blade airfoil aerodynamics model {1=steady model, 2=Beddoes-Leishman unsteady model}" - +typedef ^ AD_InputFile IntKi Wake_Mod - - - "Type of wake/induction model {0=none, 1=BEMT, 2=DBEMT, 3=FVW}" - +typedef ^ AD_InputFile IntKi BEM_Mod - - - "Type of BEM model {1=legacy NoSweepPitchTwist, 2=polar grid}" - typedef ^ AD_InputFile IntKi TwrPotent - - - "Type of tower influence on wind based on potential flow around the tower {0=none, 1=baseline potential flow, 2=potential flow with Bak correction}" - typedef ^ AD_InputFile IntKi TwrShadow - - - "Type of tower influence on wind based on downstream tower shadow {0=none, 1=Powles model, 2=Eames model}" - typedef ^ AD_InputFile LOGICAL TwrAero - - - "Calculate tower aerodynamic loads?" flag -typedef ^ AD_InputFile Logical FrozenWake - - - "Flag that tells this module it should assume a frozen wake during linearization." - typedef ^ AD_InputFile Logical CavitCheck - - - "Flag that tells us if we want to check for cavitation" - typedef ^ AD_InputFile Logical Buoyancy - - - "Include buoyancy effects?" flag typedef ^ AD_InputFile Logical CompAA - - - "Compute AeroAcoustic noise" flag @@ -192,16 +190,24 @@ typedef ^ AD_InputFile ReKi KinVisc - - - "Kinematic air viscosity" m^2/s typedef ^ AD_InputFile ReKi Patm - - - "Atmospheric pressure" Pa typedef ^ AD_InputFile ReKi Pvap - - - "Vapour pressure" Pa typedef ^ AD_InputFile ReKi SpdSound - - - "Speed of sound" m/s -typedef ^ AD_InputFile IntKi SkewMod - - - "Type of skewed-wake correction model {0=orthogonal, 1=uncoupled, 2=Pitt/Peters, 3=coupled} [unused when WakeMod=0]" - +typedef ^ AD_InputFile IntKi Skew_Mod - - - "Select skew model {0=No skew model at all, -1=Throw away non-normal component for linearization, 1=Glauert skew model}" - +typedef ^ AD_InputFile Logical SkewMomCorr - - - "Turn the skew momentum correction on or off [used only when SkewMod=1]" - +typedef ^ AD_InputFile IntKi SkewRedistr_Mod - - - "Type of skewed-wake redistribution model (switch) {0=no redistribution, 1=Glauert/Pitt/Peters, 2=Vortex Cylinder} [unsed only when SkewMod=1]" - typedef ^ AD_InputFile ReKi SkewModFactor - - - "Constant used in Pitt/Peters skewed wake model (default is 15*pi/32)" - -typedef ^ AD_InputFile LOGICAL TipLoss - - - "Use the Prandtl tip-loss model? [unused when WakeMod=0]" flag -typedef ^ AD_InputFile LOGICAL HubLoss - - - "Use the Prandtl hub-loss model? [unused when WakeMod=0]" flag -typedef ^ AD_InputFile LOGICAL TanInd - - - "Include tangential induction in BEMT calculations? [unused when WakeMod=0]" flag -typedef ^ AD_InputFile LOGICAL AIDrag - - - "Include the drag term in the axial-induction calculation? [unused when WakeMod=0]" flag -typedef ^ AD_InputFile LOGICAL TIDrag - - - "Include the drag term in the tangential-induction calculation? [unused when WakeMod=0 or TanInd=FALSE]" flag -typedef ^ AD_InputFile ReKi IndToler - - - "Convergence tolerance for BEM induction factors [unused when WakeMod=0]" - -typedef ^ AD_InputFile ReKi MaxIter - - - "Maximum number of iteration steps [unused when WakeMod=0]" - -typedef ^ AD_InputFile IntKi UAMod - - - "Unsteady Aero Model Switch (switch) {1=Baseline model (Original), 2=Gonzalez's variant (changes in Cn,Cc,Cm), 3=Minnema/Pierce variant (changes in Cc and Cm)} [used only when AFAeroMod=2]" - +typedef ^ AD_InputFile LOGICAL TipLoss - - - "Use the Prandtl tip-loss model? [unused when Wake_Mod=0]" flag +typedef ^ AD_InputFile LOGICAL HubLoss - - - "Use the Prandtl hub-loss model? [unused when Wake_Mod=0]" flag +typedef ^ AD_InputFile LOGICAL TanInd - - - "Include tangential induction in BEMT calculations? [unused when Wake_Mod=0]" flag +typedef ^ AD_InputFile LOGICAL AIDrag - - - "Include the drag term in the axial-induction calculation? [unused when Wake_Mod=0]" flag +typedef ^ AD_InputFile LOGICAL TIDrag - - - "Include the drag term in the tangential-induction calculation? [unused when Wake_Mod=0 or TanInd=FALSE]" flag +typedef ^ AD_InputFile ReKi IndToler - - - "Convergence tolerance for BEM induction factors [unused when Wake_Mod=0]" - +typedef ^ AD_InputFile ReKi MaxIter - - - "Maximum number of iteration steps [unused when Wake_Mod=0]" - +typedef ^ AD_InputFile Logical SectAvg - .False. - "Use Sector average for BEM inflow velocity calculation (flag)" - +typedef ^ ^ IntKi SA_Weighting - 1 - "Sector Average - Weighting function for sector average {1=Uniform, 2=Impulse, } within a 360/nB sector centered on the blade (switch) [used only when SectAvg=True]" - +typedef ^ ^ ReKi SA_PsiBwd - -60 - "Sector Average - Backard Azimuth (<0)" deg +typedef ^ ^ ReKi SA_PsiFwd - 60 - "Sector Average - Forward Azimuth (>0)" deg +typedef ^ ^ IntKi SA_nPerSec - 5 - "Sector average - Number of points per sectors (-) [used only when SectAvg=True]" - +typedef ^ ^ LOGICAL AoA34 - - - "Sample the angle of attack (AoA) at the 3/4 chord or the AC point {default=True} [always used]" - +typedef ^ ^ IntKi UA_Mod - - - "Unsteady Aero Model Switch (switch) {0=Quasi-steady (no UA), 2=Gonzalez's variant (changes in Cn,Cc,Cm), 3=Minnema/Pierce variant (changes in Cc and Cm)}" - typedef ^ AD_InputFile LOGICAL FLookup - - - "Flag to indicate whether a lookup for f' will be calculated (TRUE) or whether best-fit exponential equations will be used (FALSE); if FALSE S1-S4 must be provided in airfoil input files [used only when AFAeroMod=2]" flag typedef ^ AD_InputFile ReKi InCol_Alfa - - - "The column in the airfoil tables that contains the angle of attack" - typedef ^ AD_InputFile ReKi InCol_Cl - - - "The column in the airfoil tables that contains the lift coefficient" - @@ -221,8 +227,8 @@ typedef ^ AD_InputFile IntKi NTwOuts - - - "Number of tower node outputs [0 - 9] typedef ^ AD_InputFile IntKi TwOutNd {9} - - "Tower nodes whose values will be output" - typedef ^ AD_InputFile IntKi NumOuts - - - "Number of parameters in the output list (number of outputs requested)" - typedef ^ AD_InputFile CHARACTER(ChanLen) OutList {:} - - "List of user-requested output channels" - -typedef ^ AD_InputFile ReKi tau1_const - - - "time constant for DBEMT [used only when WakeMod=2 and DBEMT_Mod/=2]" s -typedef ^ AD_InputFile IntKi DBEMT_Mod - - - "Type of dynamic BEMT (DBEMT) model {1=constant tau1, 2=time-dependent tau1}" - +typedef ^ AD_InputFile ReKi tau1_const - - - "time constant for DBEMT" s +typedef ^ AD_InputFile IntKi DBEMT_Mod - - - "Type of dynamic BEMT (DBEMT) model {0=No Dynamic Wake, -1=Frozen Wake for linearization, 1=constant tau1, 2=time-dependent tau1, 3=constant tau1 with continuous formulation} (-) [used only when WakeMod=1]" - typedef ^ AD_InputFile IntKi BldNd_NumOuts - - - "Number of requested output channels per blade node (AD_AllBldNdOuts)" - typedef ^ AD_InputFile CHARACTER(ChanLen) BldNd_OutList {:} - - "List of user-requested output channels (AD_AllBldNdOuts)" - #typedef ^ AD_InputFile IntKi BldNd_BlOutNd {:} - - "The blade nodes to actually output (AD_AllBldNdOuts)" - @@ -277,6 +283,7 @@ typedef ^ RotMiscVarType AA_OutputType AA_y - - - "Outputs from the AA module" - typedef ^ RotMiscVarType AA_InputType AA_u - - - "Inputs to the AA module" - typedef ^ RotMiscVarType ReKi DisturbedInflow {:}{:}{:} - - "InflowOnBlade values modified by tower influence" m/s +typedef ^ RotMiscVarType ReKi SectAvgInflow {:}{:}{:} - - "Sector averaged - disturbed inflow to improve BEM shear calculations" m/s typedef ^ RotMiscVarType R8Ki orientationAnnulus {:}{:}{:}{:} - - "Coordinate system equivalent to BladeMotion Orientation, but without live sweep, blade-pitch, and twist angles" - typedef ^ RotMiscVarType R8Ki R_li {:}{:}{:}{:} - - "Transformation matrix from inertial system to the staggered polar coordinate system of a given section" - typedef ^ RotMiscVarType ReKi AllOuts {:} - - "An array holding the value of all of the calculated (not only selected) output channels" - @@ -379,7 +386,7 @@ typedef ^ RotParameterType Integer NumBl_Lin - - - " typedef ^ RotParameterType IntKi TwrPotent - - - "Type of tower influence on wind based on potential flow around the tower {0=none, 1=baseline potential flow, 2=potential flow with Bak correction}" - typedef ^ RotParameterType IntKi TwrShadow - - - "Type of tower influence on wind based on downstream tower shadow {0=none, 1=Powles model, 2=Eames model}" - typedef ^ RotParameterType LOGICAL TwrAero - - - "Calculate tower aerodynamic loads?" flag -typedef ^ RotParameterType Logical FrozenWake - - - "Flag that tells this module it should assume a frozen wake during linearization." - +typedef ^ RotParameterType Integer DBEMT_Mod - - - "DBEMT_Mod" - typedef ^ RotParameterType Logical CavitCheck - - - "Flag that tells us if we want to check for cavitation" - typedef ^ RotParameterType Logical Buoyancy - - - "Include buoyancy effects?" flag typedef ^ RotParameterType IntKi MHK - - - "MHK" flag @@ -393,7 +400,7 @@ typedef ^ RotParameterType ReKi Pvap - - - "Vapour pressure" P typedef ^ RotParameterType ReKi WtrDpth - - - "Water depth" m typedef ^ RotParameterType ReKi MSL2SWL - - - "Offset between still-water level and mean sea level" m typedef ^ RotParameterType IntKi AeroProjMod - 1 - "Flag to switch between different projection models" - -typedef ^ RotParameterType IntKi AeroBEM_Mod - -1 - "Flag to switch between different BEM Model" - +typedef ^ RotParameterType IntKi BEM_Mod - -1 - "Flag to switch between different BEM Model" - # parameters for output typedef ^ RotParameterType IntKi NumOuts - - - "Number of parameters in the output list (number of outputs requested)" - typedef ^ RotParameterType CHARACTER(1024) RootName - - - "RootName for writing output files" - @@ -418,12 +425,17 @@ typedef ^ ParameterType RotParameterType rotors {:} - - "Parameter types for typedef ^ ParameterType DbKi DT - - - "Time step for continuous state integration & discrete state update" seconds typedef ^ ParameterType CHARACTER(1024) RootName - - - "RootName for writing output files" - typedef ^ ParameterType AFI_ParameterType AFI {:} - - "AirfoilInfo parameters" -typedef ^ ParameterType IntKi SkewMod - - - "Type of skewed-wake correction model {0=orthogonal, 1=uncoupled, 2=Pitt/Peters, 3=coupled} [unused when WakeMod=0]" - -typedef ^ ParameterType IntKi WakeMod - - - "Type of wake/induction model {0=none, 1=BEMT, 2=DBEMT, 3=FVW}" - +typedef ^ ParameterType IntKi Skew_Mod - - - "Type of skewed-wake correction model {-1=orthogonal, 0=None, 1=Glauert} [unused when Wake_Mod=0]" - +typedef ^ ParameterType IntKi Wake_Mod - - - "Type of wake/induction model {0=none, 1=BEMT, 2=DBEMT, 3=FVW}" - typedef ^ ParameterType FVW_ParameterType FVW - - - "Parameters for FVW module" typedef ^ ParameterType LOGICAL CompAeroMaps - .FALSE. - "flag to determine if AeroDyn is computing aero maps (true) or running a normal simulation (false)" - typedef ^ ParameterType LOGICAL UA_Flag - - - "logical flag indicating whether to use UnsteadyAero" - typedef ^ ParameterType FlowFieldType *FlowField - - - "Pointer of InflowWinds flow field data type" - +typedef ^ ^ Logical SectAvg - - - "Use Sector average for BEM inflow velocity calculation" - +typedef ^ ^ IntKi SA_Weighting - - 1 "Sector Average - Weighting function for sector average {1=Uniform, 2=Impulse} within a 360/nB sector centered on the blade (switch) [used only when SectAvg=True]" - +typedef ^ ^ ReKi SA_PsiBwd - - - "Sector Average - Backard Azimuth (<0)" deg +typedef ^ ^ ReKi SA_PsiFwd - - - "Sector Average - Forward Azimuth (>0)" deg +typedef ^ ^ IntKi SA_nPerSec - - - "Sector Average - Number of points per sector (>1)" - # ..... Inputs .................................................................................................................... diff --git a/modules/aerodyn/src/AeroDyn_Types.f90 b/modules/aerodyn/src/AeroDyn_Types.f90 index 0d135968a4..407e4bca26 100644 --- a/modules/aerodyn/src/AeroDyn_Types.f90 +++ b/modules/aerodyn/src/AeroDyn_Types.f90 @@ -41,16 +41,14 @@ MODULE AeroDyn_Types INTEGER(IntKi), PUBLIC, PARAMETER :: ModelUnknown = -1 ! [-] INTEGER(IntKi), PUBLIC, PARAMETER :: WakeMod_none = 0 ! Wake model - none [-] INTEGER(IntKi), PUBLIC, PARAMETER :: WakeMod_BEMT = 1 ! Wake model - BEMT (blade elememnt momentum theory) [-] - INTEGER(IntKi), PUBLIC, PARAMETER :: WakeMod_DBEMT = 2 ! Wake model - DBEMT (dynamic elememnt momentum theory) [-] INTEGER(IntKi), PUBLIC, PARAMETER :: WakeMod_FVW = 3 ! Wake model - FVW (free vortex wake, OLAF) [-] - INTEGER(IntKi), PUBLIC, PARAMETER :: AFAeroMod_steady = 1 ! steady model [-] - INTEGER(IntKi), PUBLIC, PARAMETER :: AFAeroMod_BL_unsteady = 2 ! Beddoes-Leishman unsteady model [-] INTEGER(IntKi), PUBLIC, PARAMETER :: TwrPotent_none = 0 ! no tower potential flow [-] INTEGER(IntKi), PUBLIC, PARAMETER :: TwrPotent_baseline = 1 ! baseline tower potential flow [-] INTEGER(IntKi), PUBLIC, PARAMETER :: TwrPotent_Bak = 2 ! tower potential flow with Bak correction [-] INTEGER(IntKi), PUBLIC, PARAMETER :: TwrShadow_none = 0 ! no tower shadow [-] INTEGER(IntKi), PUBLIC, PARAMETER :: TwrShadow_Powles = 1 ! Powles tower shadow model [-] INTEGER(IntKi), PUBLIC, PARAMETER :: TwrShadow_Eames = 2 ! Eames tower shadow model [-] + INTEGER(IntKi), PUBLIC, PARAMETER :: SA_Wgt_Uniform = 1 ! Sector average weighting - Uniform [-] INTEGER(IntKi), PUBLIC, PARAMETER :: TFinAero_none = 0 ! no tail fin aero [-] INTEGER(IntKi), PUBLIC, PARAMETER :: TFinAero_polar = 1 ! polar-based tail fin aerodynamics [-] INTEGER(IntKi), PUBLIC, PARAMETER :: TFinAero_USB = 2 ! unsteady slender body tail fin aerodynamics model [-] @@ -112,7 +110,6 @@ MODULE AeroDyn_Types REAL(R8Ki) , DIMENSION(1:3) :: NacellePosition = 0.0_R8Ki !< X-Y-Z reference position of nacelle [m] REAL(R8Ki) , DIMENSION(1:3,1:3) :: NacelleOrientation = 0.0_R8Ki !< DCM reference orientation of nacelle [-] INTEGER(IntKi) :: AeroProjMod = 1 !< Flag to switch between different projection models [-] - INTEGER(IntKi) :: AeroBEM_Mod = -1 !< Flag to switch between different BEM Model [-] REAL(ReKi) :: RotSpeed = 0.0_ReKi !< Rotor speed used when AeroDyn is computing aero maps [rad/s] END TYPE RotInitInputType ! ======================= @@ -204,12 +201,11 @@ MODULE AeroDyn_Types TYPE, PUBLIC :: AD_InputFile LOGICAL :: Echo = .false. !< Echo input file to echo file [-] REAL(DbKi) :: DTAero = 0.0_R8Ki !< Time interval for aerodynamic calculations {or "default"} [s] - INTEGER(IntKi) :: WakeMod = 0_IntKi !< Type of wake/induction model {0=none, 1=BEMT, 2=DBEMT, 3=FVW} [-] - INTEGER(IntKi) :: AFAeroMod = 0_IntKi !< Type of blade airfoil aerodynamics model {1=steady model, 2=Beddoes-Leishman unsteady model} [-] + INTEGER(IntKi) :: Wake_Mod = 0_IntKi !< Type of wake/induction model {0=none, 1=BEMT, 2=DBEMT, 3=FVW} [-] + INTEGER(IntKi) :: BEM_Mod = 0_IntKi !< Type of BEM model {1=legacy NoSweepPitchTwist, 2=polar grid} [-] INTEGER(IntKi) :: TwrPotent = 0_IntKi !< Type of tower influence on wind based on potential flow around the tower {0=none, 1=baseline potential flow, 2=potential flow with Bak correction} [-] INTEGER(IntKi) :: TwrShadow = 0_IntKi !< Type of tower influence on wind based on downstream tower shadow {0=none, 1=Powles model, 2=Eames model} [-] LOGICAL :: TwrAero = .false. !< Calculate tower aerodynamic loads? [flag] - LOGICAL :: FrozenWake = .false. !< Flag that tells this module it should assume a frozen wake during linearization. [-] LOGICAL :: CavitCheck = .false. !< Flag that tells us if we want to check for cavitation [-] LOGICAL :: Buoyancy = .false. !< Include buoyancy effects? [flag] LOGICAL :: CompAA = .false. !< Compute AeroAcoustic noise [flag] @@ -220,16 +216,24 @@ MODULE AeroDyn_Types REAL(ReKi) :: Patm = 0.0_ReKi !< Atmospheric pressure [Pa] REAL(ReKi) :: Pvap = 0.0_ReKi !< Vapour pressure [Pa] REAL(ReKi) :: SpdSound = 0.0_ReKi !< Speed of sound [m/s] - INTEGER(IntKi) :: SkewMod = 0_IntKi !< Type of skewed-wake correction model {0=orthogonal, 1=uncoupled, 2=Pitt/Peters, 3=coupled} [unused when WakeMod=0] [-] + INTEGER(IntKi) :: Skew_Mod = 0_IntKi !< Select skew model {0=No skew model at all, -1=Throw away non-normal component for linearization, 1=Glauert skew model} [-] + LOGICAL :: SkewMomCorr = .false. !< Turn the skew momentum correction on or off [used only when SkewMod=1] [-] + INTEGER(IntKi) :: SkewRedistr_Mod = 0_IntKi !< Type of skewed-wake redistribution model (switch) {0=no redistribution, 1=Glauert/Pitt/Peters, 2=Vortex Cylinder} [unsed only when SkewMod=1] [-] REAL(ReKi) :: SkewModFactor = 0.0_ReKi !< Constant used in Pitt/Peters skewed wake model (default is 15*pi/32) [-] - LOGICAL :: TipLoss = .false. !< Use the Prandtl tip-loss model? [unused when WakeMod=0] [flag] - LOGICAL :: HubLoss = .false. !< Use the Prandtl hub-loss model? [unused when WakeMod=0] [flag] - LOGICAL :: TanInd = .false. !< Include tangential induction in BEMT calculations? [unused when WakeMod=0] [flag] - LOGICAL :: AIDrag = .false. !< Include the drag term in the axial-induction calculation? [unused when WakeMod=0] [flag] - LOGICAL :: TIDrag = .false. !< Include the drag term in the tangential-induction calculation? [unused when WakeMod=0 or TanInd=FALSE] [flag] - REAL(ReKi) :: IndToler = 0.0_ReKi !< Convergence tolerance for BEM induction factors [unused when WakeMod=0] [-] - REAL(ReKi) :: MaxIter = 0.0_ReKi !< Maximum number of iteration steps [unused when WakeMod=0] [-] - INTEGER(IntKi) :: UAMod = 0_IntKi !< Unsteady Aero Model Switch (switch) {1=Baseline model (Original), 2=Gonzalez's variant (changes in Cn,Cc,Cm), 3=Minnema/Pierce variant (changes in Cc and Cm)} [used only when AFAeroMod=2] [-] + LOGICAL :: TipLoss = .false. !< Use the Prandtl tip-loss model? [unused when Wake_Mod=0] [flag] + LOGICAL :: HubLoss = .false. !< Use the Prandtl hub-loss model? [unused when Wake_Mod=0] [flag] + LOGICAL :: TanInd = .false. !< Include tangential induction in BEMT calculations? [unused when Wake_Mod=0] [flag] + LOGICAL :: AIDrag = .false. !< Include the drag term in the axial-induction calculation? [unused when Wake_Mod=0] [flag] + LOGICAL :: TIDrag = .false. !< Include the drag term in the tangential-induction calculation? [unused when Wake_Mod=0 or TanInd=FALSE] [flag] + REAL(ReKi) :: IndToler = 0.0_ReKi !< Convergence tolerance for BEM induction factors [unused when Wake_Mod=0] [-] + REAL(ReKi) :: MaxIter = 0.0_ReKi !< Maximum number of iteration steps [unused when Wake_Mod=0] [-] + LOGICAL :: SectAvg = .False. !< Use Sector average for BEM inflow velocity calculation (flag) [-] + INTEGER(IntKi) :: SA_Weighting = 1 !< Sector Average - Weighting function for sector average {1=Uniform, 2=Impulse, } within a 360/nB sector centered on the blade (switch) [used only when SectAvg=True] [-] + REAL(ReKi) :: SA_PsiBwd = -60 !< Sector Average - Backard Azimuth (<0) [deg] + REAL(ReKi) :: SA_PsiFwd = 60 !< Sector Average - Forward Azimuth (>0) [deg] + INTEGER(IntKi) :: SA_nPerSec = 5 !< Sector average - Number of points per sectors (-) [used only when SectAvg=True] [-] + LOGICAL :: AoA34 = .false. !< Sample the angle of attack (AoA) at the 3/4 chord or the AC point {default=True} [always used] [-] + INTEGER(IntKi) :: UA_Mod = 0_IntKi !< Unsteady Aero Model Switch (switch) {0=Quasi-steady (no UA), 2=Gonzalez's variant (changes in Cn,Cc,Cm), 3=Minnema/Pierce variant (changes in Cc and Cm)} [-] LOGICAL :: FLookup = .false. !< Flag to indicate whether a lookup for f' will be calculated (TRUE) or whether best-fit exponential equations will be used (FALSE); if FALSE S1-S4 must be provided in airfoil input files [used only when AFAeroMod=2] [flag] REAL(ReKi) :: InCol_Alfa = 0.0_ReKi !< The column in the airfoil tables that contains the angle of attack [-] REAL(ReKi) :: InCol_Cl = 0.0_ReKi !< The column in the airfoil tables that contains the lift coefficient [-] @@ -248,8 +252,8 @@ MODULE AeroDyn_Types INTEGER(IntKi) , DIMENSION(1:9) :: TwOutNd = 0_IntKi !< Tower nodes whose values will be output [-] INTEGER(IntKi) :: NumOuts = 0_IntKi !< Number of parameters in the output list (number of outputs requested) [-] CHARACTER(ChanLen) , DIMENSION(:), ALLOCATABLE :: OutList !< List of user-requested output channels [-] - REAL(ReKi) :: tau1_const = 0.0_ReKi !< time constant for DBEMT [used only when WakeMod=2 and DBEMT_Mod/=2] [s] - INTEGER(IntKi) :: DBEMT_Mod = 0_IntKi !< Type of dynamic BEMT (DBEMT) model {1=constant tau1, 2=time-dependent tau1} [-] + REAL(ReKi) :: tau1_const = 0.0_ReKi !< time constant for DBEMT [s] + INTEGER(IntKi) :: DBEMT_Mod = 0_IntKi !< Type of dynamic BEMT (DBEMT) model {0=No Dynamic Wake, -1=Frozen Wake for linearization, 1=constant tau1, 2=time-dependent tau1, 3=constant tau1 with continuous formulation} (-) [used only when WakeMod=1] [-] INTEGER(IntKi) :: BldNd_NumOuts = 0_IntKi !< Number of requested output channels per blade node (AD_AllBldNdOuts) [-] CHARACTER(ChanLen) , DIMENSION(:), ALLOCATABLE :: BldNd_OutList !< List of user-requested output channels (AD_AllBldNdOuts) [-] CHARACTER(1024) :: BldNd_BlOutNd_Str !< String to parse for the blade nodes to actually output (AD_AllBldNdOuts) [-] @@ -317,6 +321,7 @@ MODULE AeroDyn_Types TYPE(AA_OutputType) :: AA_y !< Outputs from the AA module [-] TYPE(AA_InputType) :: AA_u !< Inputs to the AA module [-] REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: DisturbedInflow !< InflowOnBlade values modified by tower influence [m/s] + REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE :: SectAvgInflow !< Sector averaged - disturbed inflow to improve BEM shear calculations [m/s] REAL(R8Ki) , DIMENSION(:,:,:,:), ALLOCATABLE :: orientationAnnulus !< Coordinate system equivalent to BladeMotion Orientation, but without live sweep, blade-pitch, and twist angles [-] REAL(R8Ki) , DIMENSION(:,:,:,:), ALLOCATABLE :: R_li !< Transformation matrix from inertial system to the staggered polar coordinate system of a given section [-] REAL(ReKi) , DIMENSION(:), ALLOCATABLE :: AllOuts !< An array holding the value of all of the calculated (not only selected) output channels [-] @@ -418,7 +423,7 @@ MODULE AeroDyn_Types INTEGER(IntKi) :: TwrPotent = 0_IntKi !< Type of tower influence on wind based on potential flow around the tower {0=none, 1=baseline potential flow, 2=potential flow with Bak correction} [-] INTEGER(IntKi) :: TwrShadow = 0_IntKi !< Type of tower influence on wind based on downstream tower shadow {0=none, 1=Powles model, 2=Eames model} [-] LOGICAL :: TwrAero = .false. !< Calculate tower aerodynamic loads? [flag] - LOGICAL :: FrozenWake = .false. !< Flag that tells this module it should assume a frozen wake during linearization. [-] + INTEGER(IntKi) :: DBEMT_Mod = 0_IntKi !< DBEMT_Mod [-] LOGICAL :: CavitCheck = .false. !< Flag that tells us if we want to check for cavitation [-] LOGICAL :: Buoyancy = .false. !< Include buoyancy effects? [flag] INTEGER(IntKi) :: MHK = 0_IntKi !< MHK [flag] @@ -432,7 +437,7 @@ MODULE AeroDyn_Types REAL(ReKi) :: WtrDpth = 0.0_ReKi !< Water depth [m] REAL(ReKi) :: MSL2SWL = 0.0_ReKi !< Offset between still-water level and mean sea level [m] INTEGER(IntKi) :: AeroProjMod = 1 !< Flag to switch between different projection models [-] - INTEGER(IntKi) :: AeroBEM_Mod = -1 !< Flag to switch between different BEM Model [-] + INTEGER(IntKi) :: BEM_Mod = -1 !< Flag to switch between different BEM Model [-] INTEGER(IntKi) :: NumOuts = 0_IntKi !< Number of parameters in the output list (number of outputs requested) [-] CHARACTER(1024) :: RootName !< RootName for writing output files [-] TYPE(OutParmType) , DIMENSION(:), ALLOCATABLE :: OutParam !< Names and units (and other characteristics) of all requested output parameters [-] @@ -455,12 +460,17 @@ MODULE AeroDyn_Types REAL(DbKi) :: DT = 0.0_R8Ki !< Time step for continuous state integration & discrete state update [seconds] CHARACTER(1024) :: RootName !< RootName for writing output files [-] TYPE(AFI_ParameterType) , DIMENSION(:), ALLOCATABLE :: AFI !< AirfoilInfo parameters [-] - INTEGER(IntKi) :: SkewMod = 0_IntKi !< Type of skewed-wake correction model {0=orthogonal, 1=uncoupled, 2=Pitt/Peters, 3=coupled} [unused when WakeMod=0] [-] - INTEGER(IntKi) :: WakeMod = 0_IntKi !< Type of wake/induction model {0=none, 1=BEMT, 2=DBEMT, 3=FVW} [-] + INTEGER(IntKi) :: Skew_Mod = 0_IntKi !< Type of skewed-wake correction model {-1=orthogonal, 0=None, 1=Glauert} [unused when Wake_Mod=0] [-] + INTEGER(IntKi) :: Wake_Mod = 0_IntKi !< Type of wake/induction model {0=none, 1=BEMT, 2=DBEMT, 3=FVW} [-] TYPE(FVW_ParameterType) :: FVW !< Parameters for FVW module [-] LOGICAL :: CompAeroMaps = .FALSE. !< flag to determine if AeroDyn is computing aero maps (true) or running a normal simulation (false) [-] LOGICAL :: UA_Flag = .false. !< logical flag indicating whether to use UnsteadyAero [-] TYPE(FlowFieldType) , POINTER :: FlowField => NULL() !< Pointer of InflowWinds flow field data type [-] + LOGICAL :: SectAvg = .false. !< Use Sector average for BEM inflow velocity calculation [-] + INTEGER(IntKi) :: SA_Weighting = 0_IntKi !< Sector Average - Weighting function for sector average {1=Uniform, 2=Impulse} within a 360/nB sector centered on the blade (switch) [used only when SectAvg=True] [-] + REAL(ReKi) :: SA_PsiBwd = 0.0_ReKi !< Sector Average - Backard Azimuth (<0) [deg] + REAL(ReKi) :: SA_PsiFwd = 0.0_ReKi !< Sector Average - Forward Azimuth (>0) [deg] + INTEGER(IntKi) :: SA_nPerSec = 0_IntKi !< Sector Average - Number of points per sector (>1) [-] END TYPE AD_ParameterType ! ======================= ! ========= BldInputType ======= @@ -858,7 +868,6 @@ subroutine AD_CopyRotInitInputType(SrcRotInitInputTypeData, DstRotInitInputTypeD DstRotInitInputTypeData%NacellePosition = SrcRotInitInputTypeData%NacellePosition DstRotInitInputTypeData%NacelleOrientation = SrcRotInitInputTypeData%NacelleOrientation DstRotInitInputTypeData%AeroProjMod = SrcRotInitInputTypeData%AeroProjMod - DstRotInitInputTypeData%AeroBEM_Mod = SrcRotInitInputTypeData%AeroBEM_Mod DstRotInitInputTypeData%RotSpeed = SrcRotInitInputTypeData%RotSpeed end subroutine @@ -891,7 +900,6 @@ subroutine AD_PackRotInitInputType(RF, Indata) call RegPack(RF, InData%NacellePosition) call RegPack(RF, InData%NacelleOrientation) call RegPack(RF, InData%AeroProjMod) - call RegPack(RF, InData%AeroBEM_Mod) call RegPack(RF, InData%RotSpeed) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -913,7 +921,6 @@ subroutine AD_UnPackRotInitInputType(RF, OutData) call RegUnpack(RF, OutData%NacellePosition); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%NacelleOrientation); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%AeroProjMod); if (RegCheckErr(RF, RoutineName)) return - call RegUnpack(RF, OutData%AeroBEM_Mod); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%RotSpeed); if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -2012,12 +2019,11 @@ subroutine AD_CopyInputFile(SrcInputFileData, DstInputFileData, CtrlCode, ErrSta ErrMsg = '' DstInputFileData%Echo = SrcInputFileData%Echo DstInputFileData%DTAero = SrcInputFileData%DTAero - DstInputFileData%WakeMod = SrcInputFileData%WakeMod - DstInputFileData%AFAeroMod = SrcInputFileData%AFAeroMod + DstInputFileData%Wake_Mod = SrcInputFileData%Wake_Mod + DstInputFileData%BEM_Mod = SrcInputFileData%BEM_Mod DstInputFileData%TwrPotent = SrcInputFileData%TwrPotent DstInputFileData%TwrShadow = SrcInputFileData%TwrShadow DstInputFileData%TwrAero = SrcInputFileData%TwrAero - DstInputFileData%FrozenWake = SrcInputFileData%FrozenWake DstInputFileData%CavitCheck = SrcInputFileData%CavitCheck DstInputFileData%Buoyancy = SrcInputFileData%Buoyancy DstInputFileData%CompAA = SrcInputFileData%CompAA @@ -2039,7 +2045,9 @@ subroutine AD_CopyInputFile(SrcInputFileData, DstInputFileData, CtrlCode, ErrSta DstInputFileData%Patm = SrcInputFileData%Patm DstInputFileData%Pvap = SrcInputFileData%Pvap DstInputFileData%SpdSound = SrcInputFileData%SpdSound - DstInputFileData%SkewMod = SrcInputFileData%SkewMod + DstInputFileData%Skew_Mod = SrcInputFileData%Skew_Mod + DstInputFileData%SkewMomCorr = SrcInputFileData%SkewMomCorr + DstInputFileData%SkewRedistr_Mod = SrcInputFileData%SkewRedistr_Mod DstInputFileData%SkewModFactor = SrcInputFileData%SkewModFactor DstInputFileData%TipLoss = SrcInputFileData%TipLoss DstInputFileData%HubLoss = SrcInputFileData%HubLoss @@ -2048,7 +2056,13 @@ subroutine AD_CopyInputFile(SrcInputFileData, DstInputFileData, CtrlCode, ErrSta DstInputFileData%TIDrag = SrcInputFileData%TIDrag DstInputFileData%IndToler = SrcInputFileData%IndToler DstInputFileData%MaxIter = SrcInputFileData%MaxIter - DstInputFileData%UAMod = SrcInputFileData%UAMod + DstInputFileData%SectAvg = SrcInputFileData%SectAvg + DstInputFileData%SA_Weighting = SrcInputFileData%SA_Weighting + DstInputFileData%SA_PsiBwd = SrcInputFileData%SA_PsiBwd + DstInputFileData%SA_PsiFwd = SrcInputFileData%SA_PsiFwd + DstInputFileData%SA_nPerSec = SrcInputFileData%SA_nPerSec + DstInputFileData%AoA34 = SrcInputFileData%AoA34 + DstInputFileData%UA_Mod = SrcInputFileData%UA_Mod DstInputFileData%FLookup = SrcInputFileData%FLookup DstInputFileData%InCol_Alfa = SrcInputFileData%InCol_Alfa DstInputFileData%InCol_Cl = SrcInputFileData%InCol_Cl @@ -2169,12 +2183,11 @@ subroutine AD_PackInputFile(RF, Indata) if (RF%ErrStat >= AbortErrLev) return call RegPack(RF, InData%Echo) call RegPack(RF, InData%DTAero) - call RegPack(RF, InData%WakeMod) - call RegPack(RF, InData%AFAeroMod) + call RegPack(RF, InData%Wake_Mod) + call RegPack(RF, InData%BEM_Mod) call RegPack(RF, InData%TwrPotent) call RegPack(RF, InData%TwrShadow) call RegPack(RF, InData%TwrAero) - call RegPack(RF, InData%FrozenWake) call RegPack(RF, InData%CavitCheck) call RegPack(RF, InData%Buoyancy) call RegPack(RF, InData%CompAA) @@ -2185,7 +2198,9 @@ subroutine AD_PackInputFile(RF, Indata) call RegPack(RF, InData%Patm) call RegPack(RF, InData%Pvap) call RegPack(RF, InData%SpdSound) - call RegPack(RF, InData%SkewMod) + call RegPack(RF, InData%Skew_Mod) + call RegPack(RF, InData%SkewMomCorr) + call RegPack(RF, InData%SkewRedistr_Mod) call RegPack(RF, InData%SkewModFactor) call RegPack(RF, InData%TipLoss) call RegPack(RF, InData%HubLoss) @@ -2194,7 +2209,13 @@ subroutine AD_PackInputFile(RF, Indata) call RegPack(RF, InData%TIDrag) call RegPack(RF, InData%IndToler) call RegPack(RF, InData%MaxIter) - call RegPack(RF, InData%UAMod) + call RegPack(RF, InData%SectAvg) + call RegPack(RF, InData%SA_Weighting) + call RegPack(RF, InData%SA_PsiBwd) + call RegPack(RF, InData%SA_PsiFwd) + call RegPack(RF, InData%SA_nPerSec) + call RegPack(RF, InData%AoA34) + call RegPack(RF, InData%UA_Mod) call RegPack(RF, InData%FLookup) call RegPack(RF, InData%InCol_Alfa) call RegPack(RF, InData%InCol_Cl) @@ -2244,12 +2265,11 @@ subroutine AD_UnPackInputFile(RF, OutData) if (RF%ErrStat /= ErrID_None) return call RegUnpack(RF, OutData%Echo); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%DTAero); if (RegCheckErr(RF, RoutineName)) return - call RegUnpack(RF, OutData%WakeMod); if (RegCheckErr(RF, RoutineName)) return - call RegUnpack(RF, OutData%AFAeroMod); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%Wake_Mod); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%BEM_Mod); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%TwrPotent); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%TwrShadow); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%TwrAero); if (RegCheckErr(RF, RoutineName)) return - call RegUnpack(RF, OutData%FrozenWake); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%CavitCheck); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%Buoyancy); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%CompAA); if (RegCheckErr(RF, RoutineName)) return @@ -2260,7 +2280,9 @@ subroutine AD_UnPackInputFile(RF, OutData) call RegUnpack(RF, OutData%Patm); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%Pvap); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%SpdSound); if (RegCheckErr(RF, RoutineName)) return - call RegUnpack(RF, OutData%SkewMod); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%Skew_Mod); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%SkewMomCorr); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%SkewRedistr_Mod); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%SkewModFactor); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%TipLoss); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%HubLoss); if (RegCheckErr(RF, RoutineName)) return @@ -2269,7 +2291,13 @@ subroutine AD_UnPackInputFile(RF, OutData) call RegUnpack(RF, OutData%TIDrag); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%IndToler); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%MaxIter); if (RegCheckErr(RF, RoutineName)) return - call RegUnpack(RF, OutData%UAMod); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%SectAvg); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%SA_Weighting); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%SA_PsiBwd); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%SA_PsiFwd); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%SA_nPerSec); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%AoA34); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%UA_Mod); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%FLookup); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%InCol_Alfa); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%InCol_Cl); if (RegCheckErr(RF, RoutineName)) return @@ -2999,6 +3027,18 @@ subroutine AD_CopyRotMiscVarType(SrcRotMiscVarTypeData, DstRotMiscVarTypeData, C end if DstRotMiscVarTypeData%DisturbedInflow = SrcRotMiscVarTypeData%DisturbedInflow end if + if (allocated(SrcRotMiscVarTypeData%SectAvgInflow)) then + LB(1:3) = lbound(SrcRotMiscVarTypeData%SectAvgInflow, kind=B8Ki) + UB(1:3) = ubound(SrcRotMiscVarTypeData%SectAvgInflow, kind=B8Ki) + if (.not. allocated(DstRotMiscVarTypeData%SectAvgInflow)) then + allocate(DstRotMiscVarTypeData%SectAvgInflow(LB(1):UB(1),LB(2):UB(2),LB(3):UB(3)), stat=ErrStat2) + if (ErrStat2 /= 0) then + call SetErrStat(ErrID_Fatal, 'Error allocating DstRotMiscVarTypeData%SectAvgInflow.', ErrStat, ErrMsg, RoutineName) + return + end if + end if + DstRotMiscVarTypeData%SectAvgInflow = SrcRotMiscVarTypeData%SectAvgInflow + end if if (allocated(SrcRotMiscVarTypeData%orientationAnnulus)) then LB(1:4) = lbound(SrcRotMiscVarTypeData%orientationAnnulus, kind=B8Ki) UB(1:4) = ubound(SrcRotMiscVarTypeData%orientationAnnulus, kind=B8Ki) @@ -3467,6 +3507,9 @@ subroutine AD_DestroyRotMiscVarType(RotMiscVarTypeData, ErrStat, ErrMsg) if (allocated(RotMiscVarTypeData%DisturbedInflow)) then deallocate(RotMiscVarTypeData%DisturbedInflow) end if + if (allocated(RotMiscVarTypeData%SectAvgInflow)) then + deallocate(RotMiscVarTypeData%SectAvgInflow) + end if if (allocated(RotMiscVarTypeData%orientationAnnulus)) then deallocate(RotMiscVarTypeData%orientationAnnulus) end if @@ -3627,6 +3670,7 @@ subroutine AD_PackRotMiscVarType(RF, Indata) call AA_PackOutput(RF, InData%AA_y) call AA_PackInput(RF, InData%AA_u) call RegPackAlloc(RF, InData%DisturbedInflow) + call RegPackAlloc(RF, InData%SectAvgInflow) call RegPackAlloc(RF, InData%orientationAnnulus) call RegPackAlloc(RF, InData%R_li) call RegPackAlloc(RF, InData%AllOuts) @@ -3750,6 +3794,7 @@ subroutine AD_UnPackRotMiscVarType(RF, OutData) call AA_UnpackOutput(RF, OutData%AA_y) ! AA_y call AA_UnpackInput(RF, OutData%AA_u) ! AA_u call RegUnpackAlloc(RF, OutData%DisturbedInflow); if (RegCheckErr(RF, RoutineName)) return + call RegUnpackAlloc(RF, OutData%SectAvgInflow); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%orientationAnnulus); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%R_li); if (RegCheckErr(RF, RoutineName)) return call RegUnpackAlloc(RF, OutData%AllOuts); if (RegCheckErr(RF, RoutineName)) return @@ -4333,7 +4378,7 @@ subroutine AD_CopyRotParameterType(SrcRotParameterTypeData, DstRotParameterTypeD DstRotParameterTypeData%TwrPotent = SrcRotParameterTypeData%TwrPotent DstRotParameterTypeData%TwrShadow = SrcRotParameterTypeData%TwrShadow DstRotParameterTypeData%TwrAero = SrcRotParameterTypeData%TwrAero - DstRotParameterTypeData%FrozenWake = SrcRotParameterTypeData%FrozenWake + DstRotParameterTypeData%DBEMT_Mod = SrcRotParameterTypeData%DBEMT_Mod DstRotParameterTypeData%CavitCheck = SrcRotParameterTypeData%CavitCheck DstRotParameterTypeData%Buoyancy = SrcRotParameterTypeData%Buoyancy DstRotParameterTypeData%MHK = SrcRotParameterTypeData%MHK @@ -4347,7 +4392,7 @@ subroutine AD_CopyRotParameterType(SrcRotParameterTypeData, DstRotParameterTypeD DstRotParameterTypeData%WtrDpth = SrcRotParameterTypeData%WtrDpth DstRotParameterTypeData%MSL2SWL = SrcRotParameterTypeData%MSL2SWL DstRotParameterTypeData%AeroProjMod = SrcRotParameterTypeData%AeroProjMod - DstRotParameterTypeData%AeroBEM_Mod = SrcRotParameterTypeData%AeroBEM_Mod + DstRotParameterTypeData%BEM_Mod = SrcRotParameterTypeData%BEM_Mod DstRotParameterTypeData%NumOuts = SrcRotParameterTypeData%NumOuts DstRotParameterTypeData%RootName = SrcRotParameterTypeData%RootName if (allocated(SrcRotParameterTypeData%OutParam)) then @@ -4542,7 +4587,7 @@ subroutine AD_PackRotParameterType(RF, Indata) call RegPack(RF, InData%TwrPotent) call RegPack(RF, InData%TwrShadow) call RegPack(RF, InData%TwrAero) - call RegPack(RF, InData%FrozenWake) + call RegPack(RF, InData%DBEMT_Mod) call RegPack(RF, InData%CavitCheck) call RegPack(RF, InData%Buoyancy) call RegPack(RF, InData%MHK) @@ -4556,7 +4601,7 @@ subroutine AD_PackRotParameterType(RF, Indata) call RegPack(RF, InData%WtrDpth) call RegPack(RF, InData%MSL2SWL) call RegPack(RF, InData%AeroProjMod) - call RegPack(RF, InData%AeroBEM_Mod) + call RegPack(RF, InData%BEM_Mod) call RegPack(RF, InData%NumOuts) call RegPack(RF, InData%RootName) call RegPack(RF, allocated(InData%OutParam)) @@ -4633,7 +4678,7 @@ subroutine AD_UnPackRotParameterType(RF, OutData) call RegUnpack(RF, OutData%TwrPotent); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%TwrShadow); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%TwrAero); if (RegCheckErr(RF, RoutineName)) return - call RegUnpack(RF, OutData%FrozenWake); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%DBEMT_Mod); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%CavitCheck); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%Buoyancy); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%MHK); if (RegCheckErr(RF, RoutineName)) return @@ -4647,7 +4692,7 @@ subroutine AD_UnPackRotParameterType(RF, OutData) call RegUnpack(RF, OutData%WtrDpth); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%MSL2SWL); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%AeroProjMod); if (RegCheckErr(RF, RoutineName)) return - call RegUnpack(RF, OutData%AeroBEM_Mod); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%BEM_Mod); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%NumOuts); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%RootName); if (RegCheckErr(RF, RoutineName)) return if (allocated(OutData%OutParam)) deallocate(OutData%OutParam) @@ -4735,14 +4780,19 @@ subroutine AD_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) if (ErrStat >= AbortErrLev) return end do end if - DstParamData%SkewMod = SrcParamData%SkewMod - DstParamData%WakeMod = SrcParamData%WakeMod + DstParamData%Skew_Mod = SrcParamData%Skew_Mod + DstParamData%Wake_Mod = SrcParamData%Wake_Mod call FVW_CopyParam(SrcParamData%FVW, DstParamData%FVW, CtrlCode, ErrStat2, ErrMsg2) call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName) if (ErrStat >= AbortErrLev) return DstParamData%CompAeroMaps = SrcParamData%CompAeroMaps DstParamData%UA_Flag = SrcParamData%UA_Flag DstParamData%FlowField => SrcParamData%FlowField + DstParamData%SectAvg = SrcParamData%SectAvg + DstParamData%SA_Weighting = SrcParamData%SA_Weighting + DstParamData%SA_PsiBwd = SrcParamData%SA_PsiBwd + DstParamData%SA_PsiFwd = SrcParamData%SA_PsiFwd + DstParamData%SA_nPerSec = SrcParamData%SA_nPerSec end subroutine subroutine AD_DestroyParam(ParamData, ErrStat, ErrMsg) @@ -4807,8 +4857,8 @@ subroutine AD_PackParam(RF, Indata) call AFI_PackParam(RF, InData%AFI(i1)) end do end if - call RegPack(RF, InData%SkewMod) - call RegPack(RF, InData%WakeMod) + call RegPack(RF, InData%Skew_Mod) + call RegPack(RF, InData%Wake_Mod) call FVW_PackParam(RF, InData%FVW) call RegPack(RF, InData%CompAeroMaps) call RegPack(RF, InData%UA_Flag) @@ -4819,6 +4869,11 @@ subroutine AD_PackParam(RF, Indata) call IfW_FlowField_PackFlowFieldType(RF, InData%FlowField) end if end if + call RegPack(RF, InData%SectAvg) + call RegPack(RF, InData%SA_Weighting) + call RegPack(RF, InData%SA_PsiBwd) + call RegPack(RF, InData%SA_PsiFwd) + call RegPack(RF, InData%SA_nPerSec) if (RegCheckErr(RF, RoutineName)) return end subroutine @@ -4861,8 +4916,8 @@ subroutine AD_UnPackParam(RF, OutData) call AFI_UnpackParam(RF, OutData%AFI(i1)) ! AFI end do end if - call RegUnpack(RF, OutData%SkewMod); if (RegCheckErr(RF, RoutineName)) return - call RegUnpack(RF, OutData%WakeMod); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%Skew_Mod); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%Wake_Mod); if (RegCheckErr(RF, RoutineName)) return call FVW_UnpackParam(RF, OutData%FVW) ! FVW call RegUnpack(RF, OutData%CompAeroMaps); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%UA_Flag); if (RegCheckErr(RF, RoutineName)) return @@ -4884,6 +4939,11 @@ subroutine AD_UnPackParam(RF, OutData) else OutData%FlowField => null() end if + call RegUnpack(RF, OutData%SectAvg); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%SA_Weighting); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%SA_PsiBwd); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%SA_PsiFwd); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%SA_nPerSec); if (RegCheckErr(RF, RoutineName)) return end subroutine subroutine AD_CopyBldInputType(SrcBldInputTypeData, DstBldInputTypeData, CtrlCode, ErrStat, ErrMsg) diff --git a/modules/aerodyn/src/BEMT.f90 b/modules/aerodyn/src/BEMT.f90 index b99c3a6658..6be1f378ad 100644 --- a/modules/aerodyn/src/BEMT.f90 +++ b/modules/aerodyn/src/BEMT.f90 @@ -90,6 +90,21 @@ real(ReKi) function ComputePhiWithInduction( Vx, Vy, a, aprime, cantAngle, xVelC end function ComputePhiWithInduction +subroutine ComputePhiFromInductions(u, p, phi, axInduction, tanInduction) + type(BEMT_InputType), intent(in ) :: u + type(BEMT_ParameterType), intent(in ) :: p + real(ReKi), intent(inout) :: phi(:,:) + real(ReKi), intent(in ) :: axInduction(:,:) + real(ReKi), intent(in ) :: tanInduction(:,:) + integer(IntKi) :: i, j + do j = 1,p%numBlades ! Loop through all blades + do i = 1,p%numBladeNodes ! Loop through the blade nodes / elements + phi(i,j) = ComputePhiWithInduction( u%Vx(i,j), u%Vy(i,j), axInduction(i,j), tanInduction(i,j), u%cantAngle(i,j), u%xVelCorr(i,j) ) + enddo ! I - Blade nodes / elements + enddo ! J - All blades +end subroutine ComputePhiFromInductions + + !---------------------------------------------------------------------------------------------------------------------------------- subroutine BEMT_Set_UA_InitData( InitInp, interval, Init_UA_Data, errStat, errMsg ) ! This routine is called from BEMT_Init. @@ -175,10 +190,9 @@ subroutine BEMT_SetParameters( InitInp, p, errStat, errMsg ) p%numBlades = InitInp%numBlades p%UA_Flag = InitInp%UA_Flag p%DBEMT_Mod = InitInp%DBEMT_Mod - p%MomentumCorr = InitInp%MomentumCorr p%BEM_Mod = InitInp%BEM_Mod !call WrScr('>>>> BEM_Mod '//trim(num2lstr(p%BEM_Mod))) - if ((p%BEM_Mod/=BEMMod_2D .and. p%BEM_Mod/=BEMMod_3D )) then + if (.not.(ANY( p%BEM_Mod == (/BEMMod_2D, BEMMod_3D/)))) then call SetErrStat( ErrID_Fatal, 'BEM_Mod needs to be 0 or 2 for now', errStat, errMsg, RoutineName ) return endif @@ -238,6 +252,13 @@ subroutine BEMT_SetParameters( InitInp, p, errStat, errMsg ) p%airDens = InitInp%airDens p%kinVisc = InitInp%kinVisc p%skewWakeMod = InitInp%skewWakeMod + if (p%skewWakeMod==Skew_Mod_Active) then + p%SkewRedistrMod = InitInp%SkewRedistrMod + p%MomentumCorr = InitInp%MomentumCorr + else + p%SkewRedistrMod = SkewRedistrMod_None + p%MomentumCorr = .false. + endif p%yawCorrFactor = InitInp%yawCorrFactor p%useTipLoss = InitInp%useTipLoss p%useHubLoss = InitInp%useHubLoss @@ -584,7 +605,8 @@ subroutine BEMT_Init( InitInp, u, p, x, xd, z, OtherState, AFInfo, y, misc, Inte call BEMT_InitOtherStates( OtherState, p, errStat, errMsg ) ! initialize the other states if (errStat >= AbortErrLev) return - if ( p%DBEMT_Mod /= DBEMT_none ) then + InitInp_DBEMT%DBEMT_Mod = p%DBEMT_Mod + if ( p%DBEMT_Mod > DBEMT_none ) then InitInp_DBEMT%DBEMT_Mod = p%DBEMT_Mod InitInp_DBEMT%numBlades = p%numBlades InitInp_DBEMT%numNodes = p%numBladeNodes @@ -723,7 +745,7 @@ subroutine BEMT_ReInit(p,x,xd,z,OtherState,misc,ErrStat,ErrMsg) if (p%UseInduction) then OtherState%ValidPhi = .true. - if (p%DBEMT_Mod /= DBEMT_none ) then + if (p%DBEMT_Mod > DBEMT_none ) then call DBEMT_ReInit(p%DBEMT, x%DBEMT, OtherState%DBEMT, misc%DBEMT) end if @@ -1269,7 +1291,8 @@ subroutine BEMT_CalcOutput( t, u, p, x, xd, z, OtherState, AFInfo, y, m, errStat ! calculate inductions using BEMT, applying the DBEMT, and/or skewed wake corrections as applicable: ! NOTE that we don't use the DBEMT inputs when calling its CalcOutput routine, so we'll skip calculating them here !............................................ - call BEMT_CalcOutput_Inductions( InputIndex, t, .false., .true., y%phi, u, p, x, xd, z, OtherState, AFInfo, y%axInduction, y%tanInduction, y%chi, m, errStat, errMsg ) + call BEMT_CalcOutput_Inductions( InputIndex, t, .false., .true., y%phi, u, p, x, xd, z, OtherState, AFInfo, y%axInduction, y%tanInduction, y%chi, m, errStat, errMsg,& + y%axInduction_qs, y%tanInduction_qs, y%k, y%k_p, y%F) !............................................ ! update phi if necessary (consistent with inductions) and calculate inputs to UA (EVEN if UA isn't used, because we use the inputs later): @@ -1389,7 +1412,7 @@ subroutine BEMT_InitStates(t, u, p, x, xd, z, OtherState, m, AFInfo, ErrStat, Er m%phi = z%phi call BEMT_CalcOutput_Inductions( InputIndex, t, CalculateDBEMTInputs, ApplyCorrections, m%phi, u, p, x, xd, z, OtherState, AFInfo, m%axInduction, m%tanInduction, m%chi, m, errStat, errMsg ) - if (p%DBEMT_Mod /= DBEMT_none) then + if (p%DBEMT_Mod > DBEMT_none) then call DBEMT_InitStates_AllNodes( m%u_DBEMT(InputIndex), p%DBEMT, x%DBEMT, OtherState%DBEMT ) end if @@ -1493,7 +1516,7 @@ subroutine BEMT_CalcOutput_Inductions( InputIndex, t, CalculateDBEMTInputs, Appl !............................................ ! apply DBEMT correction to axInduction and tanInduction: !............................................ - if (p%DBEMT_Mod /= DBEMT_none) then + if (p%DBEMT_Mod > DBEMT_none) then ! If we are using DBEMT, then we will obtain the time-filtered versions of axInduction(i,j), tanInduction(i,j) ! Note that the outputs of DBEMT are the state variables x%vind, so we don't NEED to set the inputs except on initialization step (when we output the inputs instead of the states) @@ -1576,7 +1599,7 @@ subroutine ApplySkewedWakeCorrection_AllNodes(p, u, m, x, phi, OtherState, axInd !............................................ ! Apply skewed wake correction to the axial induction (y%axInduction) !............................................ - if ( p%skewWakeMod == SkewMod_PittPeters ) then + if ( p%skewWakeMod == Skew_Mod_Active ) then if (p%BEM_Mod==BEMMod_2D) then ! do nothing else @@ -1588,7 +1611,7 @@ subroutine ApplySkewedWakeCorrection_AllNodes(p, u, m, x, phi, OtherState, axInd do i = 1,p%numBladeNodes ! Loop through the blade nodes / elements if ( .not. p%FixedInductions(i,j) ) then F = getHubTipLossCorrection(p%BEM_Mod, p%useHubLoss, p%useTipLoss, p%hubLossConst(i,j), p%tipLossConst(i,j), phi(i,j), u%cantAngle(i,j) ) - call ApplySkewedWakeCorrection( p%BEM_Mod, p%skewWakeMod, p%yawCorrFactor, F, u%psi_s(j), u%psiSkewOffset, u%chi0, u%rlocal(i,j)/m%Rtip(j), axInduction(i,j), chi(i,j), m%FirstWarn_Skew ) + call ApplySkewedWakeCorrection( p%BEM_Mod, p%SkewRedistrMod, p%yawCorrFactor, F, u%psi_s(j), u%psiSkewOffset, u%chi0, u%rlocal(i,j)/m%Rtip(j), axInduction(i,j), chi(i,j), m%FirstWarn_Skew ) end if ! .not. p%FixedInductions (special case for tip and/or hub loss) enddo ! I - Blade nodes / elements enddo ! J - All blades @@ -1676,7 +1699,7 @@ subroutine BEMT_CalcContStateDeriv( t, u, p, x, xd, z, OtherState, m, dxdt, AFIn !............................................................................................................................... ! compute derivatives for DBEMT continuous states: !............................................................................................................................... - if (p%DBEMT_Mod /= DBEMT_none) then + if (p%DBEMT_Mod > DBEMT_none) then if (.not. allocated(dxdt%DBEMT%element)) then call DBEMT_CopyContState( x%DBEMT, dxdt%DBEMT, MESH_UPDATECOPY, ErrStat2, ErrMsg2 ) call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName) diff --git a/modules/aerodyn/src/BEMTUncoupled.f90 b/modules/aerodyn/src/BEMTUncoupled.f90 index bd825bc739..ca71bc4dc1 100644 --- a/modules/aerodyn/src/BEMTUncoupled.f90 +++ b/modules/aerodyn/src/BEMTUncoupled.f90 @@ -59,6 +59,8 @@ module BEMTUnCoupled public :: GetEulerAnglesFromOrientation public :: VelocityIsZero + + public :: BEMTU_Test_ACT_Relationship contains !.................................................................................................................................. @@ -227,6 +229,20 @@ subroutine computeAirfoilOperatingAOA( BEM_Mod, phi, theta, cantAngle, toeAngle, end subroutine computeAirfoilOperatingAOA +! --- +!> Angle of attack in the airfoil reference frame +real(ReKi) function computeAirfoilAOA(Vrel_a) result(AoA) + real(ReKi), intent(in ) :: Vrel_a(3) + real(ReKi) :: numer, denom, ratio + ! Determine angle of attack as angle between airfoil chordline (afAxialVec) and inflow + numer = Vrel_a(2) + denom = sqrt(Vrel_a(1)**2 + Vrel_a(2)**2) + ratio = numer / denom + AoA = acos( max( min( ratio, 1.0_ReKi ), -1.0_ReKi ) ) + AoA = sign( AoA, Vrel_a(1) ) +end function computeAirfoilAOA + + !.................................................................................................................................. !> Transform the aerodynamic coefficients (Cl,Cd,Cm) (directed based on Vrel_xy_a ) !! from the airfoil coordinate system (a) to the without sweep pitch coordinate system (w) @@ -433,10 +449,10 @@ real(ReKi) function BEMTU_InductionWithResidual(p, u, i, j, phi, AFInfo, IsValid end function BEMTU_InductionWithResidual !----------------------------------------------------------------------------------------- -subroutine ApplySkewedWakeCorrection(BEM_Mod, SkewMod, yawCorrFactor, F, azimuth, azimuthOffset, chi0, tipRatio, a, chi, FirstWarn ) +subroutine ApplySkewedWakeCorrection(BEM_Mod, SkewRedistrMod, yawCorrFactor, F, azimuth, azimuthOffset, chi0, tipRatio, a, chi, FirstWarn ) integer(IntKi), intent(in ) :: BEM_Mod - integer(IntKi), intent(in ) :: SkewMod + integer(IntKi), intent(in ) :: SkewRedistrMod real(ReKi), intent(in ) :: yawCorrFactor ! set to 15*pi/32 previously; now allowed to be input (to better match data) real(ReKi), intent(in ) :: F ! tip/hub loss factor real(ReKi), intent(in ) :: azimuth @@ -450,7 +466,10 @@ subroutine ApplySkewedWakeCorrection(BEM_Mod, SkewMod, yawCorrFactor, F, azimuth ! Local variables real(ReKi) :: yawCorr real(ReKi) :: yawCorr_tan ! magnitude of the tan(chi/2) correction term (with possible limits) - + + if (SkewRedistrMod==SkewRedistrMod_None) then + return + endif ! Skewed wake correction if(BEM_Mod==BEMMod_2D) then @@ -465,7 +484,7 @@ subroutine ApplySkewedWakeCorrection(BEM_Mod, SkewMod, yawCorrFactor, F, azimuth if (FirstWarn) then call WrScr( 'Warning: SkewedWakeCorrection encountered a large value of chi ('//trim(num2lstr(chi*R2D))// & - ' deg), so the yaw correction will be limited. This warning will not be repeated though the condition may persist. See the AD15 chi output channels, and'// & + ' deg), so the yaw correction will be limited. This warning will not be repeated though the condition may persist. See the AD chi output channels, and'// & ' consider turning off the Pitt/Peters skew model (set SkewMod=1) if this condition persists.'//NewLine) FirstWarn = .false. end if @@ -476,7 +495,14 @@ subroutine ApplySkewedWakeCorrection(BEM_Mod, SkewMod, yawCorrFactor, F, azimuth end if !bjj: modified 22-Sep-2015: RRD recommends 32 instead of 64 in the denominator (like AD14) - yawCorr = ( yawCorrFactor * yawCorr_tan * (tipRatio) * sin(azimuth) ) ! bjj: note that when chi gets close to +/-pi this blows up + ! TODO TODO TODO + if(BEM_Mod==BEMMod_2D) then + ! ADLEG: + yawCorr = ( yawCorrFactor * yawCorr_tan * (tipRatio) * sin(azimuth) ) ! bjj: note that when chi gets close to +/-pi this blows up + else + ! ADENV: + yawCorr = ( yawCorrFactor * F * yawCorr_tan * (tipRatio) * cos(azimuth-azimuthOffset) ) ! bjj: note that when chi gets close to +/-pi this blows up + endif a = a * (1.0 + yawCorr) @@ -1216,4 +1242,52 @@ FUNCTION GetEulerAnglesFromOrientation(EulerDCM,orientation) RESULT(theta) end function !----------------------------------------------------------------------------------------- + + +!> Simple test for a-Ct relationship. +subroutine BEMTU_Test_ACT_Relationship() + real(R8Ki) :: chi0 + real(R8Ki) :: delta_chi + real(ReKi) :: F + logical :: skewConvention + integer :: i + integer :: iUnit + real(R8Ki) :: c2, c1, c0 ! Empirical CT = c2*a^2 + c1*a + c0 for a > a0 + ! Get Coefficients for Empirical CT + iUnit = 123 + + ! --- No Momentum Corr, F=1 + F=1; skewConvention=.False. + call parametricStudy('ACTCoeffs_F10_NoCo.csv') + ! --- No Momentum Corr, F=0.5 + F=0.5; skewConvention=.False. + call parametricStudy('ACTCoeffs_F05_NoCo.csv') + ! --- Momentum Corr, F=1 + F=1; skewConvention=.True. + call parametricStudy('ACTCoeffs_F10_Corr.csv') + ! --- Momentum Corr, F=0.5 + F=0.5; skewConvention=.True. + call parametricStudy('ACTCoeffs_F05_Corr.csv') + + STOP + +contains + subroutine parametricStudy(filename) + character(len=*) :: filename + chi0=-50 * D2R + open(unit=iUnit, file=filename) + write(iUnit, '(5(A15))') 'chi0', 'c0', 'c1', 'c2', 'F' + do i=1,21 + call getEmpiricalCoefficients(chi0 ,F , c0, c1, c2, skewConvention) + write(iUnit,'(5(F15.5))') chi0*R2D, c0, c1, c2, F + chi0 = chi0 + 5*D2R + enddo + close(iUnit) + end subroutine + + + +end subroutine BEMTU_Test_ACT_Relationship + + end module BEMTUncoupled diff --git a/modules/aerodyn/src/BEMT_Registry.txt b/modules/aerodyn/src/BEMT_Registry.txt index 893fc7fdda..7119b9261a 100644 --- a/modules/aerodyn/src/BEMT_Registry.txt +++ b/modules/aerodyn/src/BEMT_Registry.txt @@ -16,13 +16,16 @@ usefrom AirfoilInfo_Registry.txt usefrom UnsteadyAero_Registry.txt usefrom DBEMT_Registry.txt -param BEMT/BEMT - INTEGER SkewMod_Orthogonal - 0 - "Inflow orthogonal to rotor [-]" - -param BEMT/BEMT - INTEGER SkewMod_Uncoupled - 1 - "Uncoupled (no correction)" - -param BEMT/BEMT - INTEGER SkewMod_PittPeters - 2 - "Pitt/Peters" - -param BEMT/BEMT - INTEGER SkewMod_Coupled - 3 - "Coupled" - -param BEMT/BEMT - INTEGER SkewMod_PittPeters_Cont - 4 - "Pitt/Peters continuous formulation" - +param BEMT/BEMT - INTEGER Skew_Mod_Orthogonal - -1 - "Inflow orthogonal to rotor [-]" - +param BEMT/BEMT - INTEGER Skew_Mod_None - 0 - "No skew model" - +param BEMT/BEMT - INTEGER Skew_Mod_Active - 1 - "Skew model active" - +param BEMT/BEMT - INTEGER Skew_Mod_PittPeters_Cont - 4 - "Pitt/Peters continuous formulation" - -param BEMT/BEMT - INTEGER BEMMod_2D - 0 - "2D BEM assuming Cx, Cy, phi, L, D are in the same plane" - +param BEMT/BEMT - INTEGER SkewRedistrMod_None - 0 - "No redistribution" - +param BEMT/BEMT - INTEGER SkewRedistrMod_PittPeters - 1 - "Pitt/Peters/Glauert redistribution" - +#param BEMT/BEMT - INTEGER SkewRedistrMod_VCyl - 2 - "Vortex cylinder redistribution" - + +param BEMT/BEMT - INTEGER BEMMod_2D - 1 - "2D BEM assuming Cx, Cy, phi, L, D are in the same plane" - param BEMT/BEMT - INTEGER BEMMod_3D - 2 - "3D BEM assuming a momentum balance system, and an airfoil system" - # @@ -35,7 +38,8 @@ typedef BEMT/BEMT InitInputType ReKi typedef ^ ^ INTEGER numBlades - - - "Number of blades" - typedef ^ ^ ReKi airDens - - - "Air density" kg/m^3 typedef ^ ^ ReKi kinVisc - - - "Kinematic air viscosity" m^2/s -typedef ^ ^ INTEGER skewWakeMod - - - "Type of skewed-wake correction model [switch] {1=uncoupled, 2=Pitt/Peters, 3=coupled}" - +typedef ^ ^ INTEGER skewWakeMod - - - "Type of skewed-wake model [switch] {0=None, 1=Glauert}" - +typedef ^ ^ INTEGER skewRedistrMod - - - "Type of skewed-wake redistribution model (switch) {0=no redistribution, 1=Glauert/Pitt/Peters, 2=Vortex Cylinder} [unsed only when SkewMod=1]" - typedef ^ ^ ReKi aTol - - - "Tolerance for the induction solution" - typedef ^ ^ LOGICAL useTipLoss - - - "Use the Prandtl tip-loss model? [flag]" - typedef ^ ^ LOGICAL useHubLoss - - - "Use the Prandtl hub-loss model? [flag]" - @@ -139,7 +143,8 @@ typedef ^ ^ ReKi typedef ^ ^ INTEGER numBlades - - - "Number of blades" - typedef ^ ^ ReKi airDens - - - "Air density" kg/m^3 typedef ^ ^ ReKi kinVisc - - - "Kinematic air viscosity" m^2/s -typedef ^ ^ INTEGER skewWakeMod - - - "Type of skewed-wake correction model [switch] {1=uncoupled, 2=Pitt/Peters, 3=coupled}" - +typedef ^ ^ INTEGER skewWakeMod - - - "Type of skewed-wake correction model [switch] {0=None, 1=Glauert/Pitt/Peters}" - +typedef ^ ^ INTEGER skewRedistrMod - - - "Type of skewed-wake redistribution model (switch) {0=no redistribution, 1=Glauert/Pitt/Peters, 2=Vortex Cylinder} [unsed only when SkewMod=1]" - typedef ^ ^ ReKi aTol - - - "Tolerance for the induction solution" - typedef ^ ^ LOGICAL useTipLoss - - - "Use the Prandtl tip-loss model? [flag]" - typedef ^ ^ LOGICAL useHubLoss - - - "Use the Prandtl hub-loss model? [flag]" - diff --git a/modules/aerodyn/src/BEMT_Types.f90 b/modules/aerodyn/src/BEMT_Types.f90 index 01b8ef02a6..b4ae2d590a 100644 --- a/modules/aerodyn/src/BEMT_Types.f90 +++ b/modules/aerodyn/src/BEMT_Types.f90 @@ -36,12 +36,13 @@ MODULE BEMT_Types USE DBEMT_Types USE NWTC_Library IMPLICIT NONE - INTEGER(IntKi), PUBLIC, PARAMETER :: SkewMod_Orthogonal = 0 ! Inflow orthogonal to rotor [-] [-] - INTEGER(IntKi), PUBLIC, PARAMETER :: SkewMod_Uncoupled = 1 ! Uncoupled (no correction) [-] - INTEGER(IntKi), PUBLIC, PARAMETER :: SkewMod_PittPeters = 2 ! Pitt/Peters [-] - INTEGER(IntKi), PUBLIC, PARAMETER :: SkewMod_Coupled = 3 ! Coupled [-] - INTEGER(IntKi), PUBLIC, PARAMETER :: SkewMod_PittPeters_Cont = 4 ! Pitt/Peters continuous formulation [-] - INTEGER(IntKi), PUBLIC, PARAMETER :: BEMMod_2D = 0 ! 2D BEM assuming Cx, Cy, phi, L, D are in the same plane [-] + INTEGER(IntKi), PUBLIC, PARAMETER :: Skew_Mod_Orthogonal = -1 ! Inflow orthogonal to rotor [-] [-] + INTEGER(IntKi), PUBLIC, PARAMETER :: Skew_Mod_None = 0 ! No skew model [-] + INTEGER(IntKi), PUBLIC, PARAMETER :: Skew_Mod_Active = 1 ! Skew model active [-] + INTEGER(IntKi), PUBLIC, PARAMETER :: Skew_Mod_PittPeters_Cont = 4 ! Pitt/Peters continuous formulation [-] + INTEGER(IntKi), PUBLIC, PARAMETER :: SkewRedistrMod_None = 0 ! No redistribution [-] + INTEGER(IntKi), PUBLIC, PARAMETER :: SkewRedistrMod_PittPeters = 1 ! Pitt/Peters/Glauert redistribution [-] + INTEGER(IntKi), PUBLIC, PARAMETER :: BEMMod_2D = 1 ! 2D BEM assuming Cx, Cy, phi, L, D are in the same plane [-] INTEGER(IntKi), PUBLIC, PARAMETER :: BEMMod_3D = 2 ! 3D BEM assuming a momentum balance system, and an airfoil system [-] ! ========= BEMT_InitInputType ======= TYPE, PUBLIC :: BEMT_InitInputType @@ -49,7 +50,8 @@ MODULE BEMT_Types INTEGER(IntKi) :: numBlades = 0_IntKi !< Number of blades [-] REAL(ReKi) :: airDens = 0.0_ReKi !< Air density [kg/m^3] REAL(ReKi) :: kinVisc = 0.0_ReKi !< Kinematic air viscosity [m^2/s] - INTEGER(IntKi) :: skewWakeMod = 0_IntKi !< Type of skewed-wake correction model [switch] {1=uncoupled, 2=Pitt/Peters, 3=coupled} [-] + INTEGER(IntKi) :: skewWakeMod = 0_IntKi !< Type of skewed-wake model [switch] {0=None, 1=Glauert} [-] + INTEGER(IntKi) :: skewRedistrMod = 0_IntKi !< Type of skewed-wake redistribution model (switch) {0=no redistribution, 1=Glauert/Pitt/Peters, 2=Vortex Cylinder} [unsed only when SkewMod=1] [-] REAL(ReKi) :: aTol = 0.0_ReKi !< Tolerance for the induction solution [-] LOGICAL :: useTipLoss = .false. !< Use the Prandtl tip-loss model? [flag] [-] LOGICAL :: useHubLoss = .false. !< Use the Prandtl hub-loss model? [flag] [-] @@ -150,7 +152,8 @@ MODULE BEMT_Types INTEGER(IntKi) :: numBlades = 0_IntKi !< Number of blades [-] REAL(ReKi) :: airDens = 0.0_ReKi !< Air density [kg/m^3] REAL(ReKi) :: kinVisc = 0.0_ReKi !< Kinematic air viscosity [m^2/s] - INTEGER(IntKi) :: skewWakeMod = 0_IntKi !< Type of skewed-wake correction model [switch] {1=uncoupled, 2=Pitt/Peters, 3=coupled} [-] + INTEGER(IntKi) :: skewWakeMod = 0_IntKi !< Type of skewed-wake correction model [switch] {0=None, 1=Glauert/Pitt/Peters} [-] + INTEGER(IntKi) :: skewRedistrMod = 0_IntKi !< Type of skewed-wake redistribution model (switch) {0=no redistribution, 1=Glauert/Pitt/Peters, 2=Vortex Cylinder} [unsed only when SkewMod=1] [-] REAL(ReKi) :: aTol = 0.0_ReKi !< Tolerance for the induction solution [-] LOGICAL :: useTipLoss = .false. !< Use the Prandtl tip-loss model? [flag] [-] LOGICAL :: useHubLoss = .false. !< Use the Prandtl hub-loss model? [flag] [-] @@ -256,6 +259,7 @@ subroutine BEMT_CopyInitInput(SrcInitInputData, DstInitInputData, CtrlCode, ErrS DstInitInputData%airDens = SrcInitInputData%airDens DstInitInputData%kinVisc = SrcInitInputData%kinVisc DstInitInputData%skewWakeMod = SrcInitInputData%skewWakeMod + DstInitInputData%skewRedistrMod = SrcInitInputData%skewRedistrMod DstInitInputData%aTol = SrcInitInputData%aTol DstInitInputData%useTipLoss = SrcInitInputData%useTipLoss DstInitInputData%useHubLoss = SrcInitInputData%useHubLoss @@ -421,6 +425,7 @@ subroutine BEMT_PackInitInput(RF, Indata) call RegPack(RF, InData%airDens) call RegPack(RF, InData%kinVisc) call RegPack(RF, InData%skewWakeMod) + call RegPack(RF, InData%skewRedistrMod) call RegPack(RF, InData%aTol) call RegPack(RF, InData%useTipLoss) call RegPack(RF, InData%useHubLoss) @@ -466,6 +471,7 @@ subroutine BEMT_UnPackInitInput(RF, OutData) call RegUnpack(RF, OutData%airDens); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%kinVisc); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%skewWakeMod); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%skewRedistrMod); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%aTol); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%useTipLoss); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%useHubLoss); if (RegCheckErr(RF, RoutineName)) return @@ -1214,6 +1220,7 @@ subroutine BEMT_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg) DstParamData%airDens = SrcParamData%airDens DstParamData%kinVisc = SrcParamData%kinVisc DstParamData%skewWakeMod = SrcParamData%skewWakeMod + DstParamData%skewRedistrMod = SrcParamData%skewRedistrMod DstParamData%aTol = SrcParamData%aTol DstParamData%useTipLoss = SrcParamData%useTipLoss DstParamData%useHubLoss = SrcParamData%useHubLoss @@ -1358,6 +1365,7 @@ subroutine BEMT_PackParam(RF, Indata) call RegPack(RF, InData%airDens) call RegPack(RF, InData%kinVisc) call RegPack(RF, InData%skewWakeMod) + call RegPack(RF, InData%skewRedistrMod) call RegPack(RF, InData%aTol) call RegPack(RF, InData%useTipLoss) call RegPack(RF, InData%useHubLoss) @@ -1400,6 +1408,7 @@ subroutine BEMT_UnPackParam(RF, OutData) call RegUnpack(RF, OutData%airDens); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%kinVisc); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%skewWakeMod); if (RegCheckErr(RF, RoutineName)) return + call RegUnpack(RF, OutData%skewRedistrMod); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%aTol); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%useTipLoss); if (RegCheckErr(RF, RoutineName)) return call RegUnpack(RF, OutData%useHubLoss); if (RegCheckErr(RF, RoutineName)) return diff --git a/modules/aerodyn/src/DBEMT_Registry.txt b/modules/aerodyn/src/DBEMT_Registry.txt index 21a726a068..051eb52655 100644 --- a/modules/aerodyn/src/DBEMT_Registry.txt +++ b/modules/aerodyn/src/DBEMT_Registry.txt @@ -12,6 +12,7 @@ # ...... Include files (definitions from NWTC Library) ............................................................................ include Registry_NWTC_Library.txt +param DBEMT/DBEMT - INTEGER DBEMT_frozen - -1 - "use frozen-wake for linearization (not DBEMT)" - param DBEMT/DBEMT - INTEGER DBEMT_none - 0 - "use BEMT instead (not DBEMT)" - param DBEMT/DBEMT - INTEGER DBEMT_tauConst - 1 - "use constant tau1" - param DBEMT/DBEMT - INTEGER DBEMT_tauVaries - 2 - "use time-dependent tau1" - diff --git a/modules/aerodyn/src/DBEMT_Types.f90 b/modules/aerodyn/src/DBEMT_Types.f90 index 17d9640ed6..9721c81013 100644 --- a/modules/aerodyn/src/DBEMT_Types.f90 +++ b/modules/aerodyn/src/DBEMT_Types.f90 @@ -33,6 +33,7 @@ MODULE DBEMT_Types !--------------------------------------------------------------------------------------------------------------------------------- USE NWTC_Library IMPLICIT NONE + INTEGER(IntKi), PUBLIC, PARAMETER :: DBEMT_frozen = -1 ! use frozen-wake for linearization (not DBEMT) [-] INTEGER(IntKi), PUBLIC, PARAMETER :: DBEMT_none = 0 ! use BEMT instead (not DBEMT) [-] INTEGER(IntKi), PUBLIC, PARAMETER :: DBEMT_tauConst = 1 ! use constant tau1 [-] INTEGER(IntKi), PUBLIC, PARAMETER :: DBEMT_tauVaries = 2 ! use time-dependent tau1 [-] diff --git a/modules/aerodyn/src/FVW.f90 b/modules/aerodyn/src/FVW.f90 index e501b2d9c7..1330a461fc 100644 --- a/modules/aerodyn/src/FVW.f90 +++ b/modules/aerodyn/src/FVW.f90 @@ -481,7 +481,7 @@ subroutine FVW_FinalWrite(u, p, x, z, m, ErrStat, ErrMsg) ErrMsg = "" ! Place any last minute operations or calculations here: if (p%WrVTK>0 .and. m%VTKstep>> FINAL WRITE' + call WrScr('OLAF: writting final VTK outputs') t=-1.0_ReKi if (p%WrVTK==1) then if (m%VTKstep This subroutine parses a comment line - SUBROUTINE ParseCom ( FileInfo, LineNum, Var, ErrStat, ErrMsg, UnEc ) + SUBROUTINE ParseCom ( FileInfo, LineNum, Var, ErrStat, ErrMsg, UnEc, IsLegalComment ) ! Arguments declarations. INTEGER(IntKi), INTENT(OUT) :: ErrStat !< The error status. @@ -2968,6 +2968,7 @@ SUBROUTINE ParseCom ( FileInfo, LineNum, Var, ErrStat, ErrMsg, UnEc ) CHARACTER(*), INTENT(OUT) :: Var !< The variable to receive the comment CHARACTER(*), INTENT(OUT) :: ErrMsg !< The error message, if ErrStat /= 0. TYPE (FileInfoType), INTENT(IN) :: FileInfo !< The derived type for holding the file information. + LOGICAL, OPTIONAL, INTENT(INOUT) :: IsLegalComment !< True if the comment is a "legal" comment line starting with '---' or '==='. NOTE: We have too many options, we need to be more strict!!!! CHARACTER(*), PARAMETER :: RoutineName = 'ParseCom' ErrStat=ErrID_None @@ -2995,6 +2996,21 @@ SUBROUTINE ParseCom ( FileInfo, LineNum, Var, ErrStat, ErrMsg, UnEc ) END IF LineNum = LineNum + 1 + IF (PRESENT(IsLegalComment) ) then + if (len(Var)<=3) then + IsLegalComment=.False. + else + ! Here, we are talking about comments in the input file that are "expected to be there" + IsLegalComment = (Var(1:3)=='---') .or. (Var(1:3)=='===') + endif + if (.not.IsLegalComment) then + call SetErrStat(ErrID_Fatal, NewLine//' >> A fatal error occurred when parsing data.'//NewLine// & + ' >> The comment line did not start with `---` or `===`. LineNum='// & + trim(num2lstr(LineNum))//'; NumLines='//trim(num2lstr(size(FileInfo%Lines))) & + , ErrStat, ErrMsg, RoutineName ) + endif + END IF + END SUBROUTINE ParseCom !======================================================================= diff --git a/reg_tests/r-test b/reg_tests/r-test index 5a064f10ac..8e5488c57e 160000 --- a/reg_tests/r-test +++ b/reg_tests/r-test @@ -1 +1 @@ -Subproject commit 5a064f10ac66f0808f4d0ffb0973c1c15b44c6c3 +Subproject commit 8e5488c57e06fb191cc28cb5a32e7b5cfda5a5ea