Merge remote-tracking branch 'OpenFAST/dev' into f/ua-elast-dev

docs/conf.py

@@ -72,6 +72,7 @@ def runDoxygen(sourcfile, doxyfileIn, doxyfileOut):
     'source/user/aerodyn-aeroacoustics/references.bib',
     'source/user/aerodyn-olaf/bibliography.bib',
     'source/user/aerodyn/bibliography.bib',
+    'source/user/elastodyn/bibliography.bib',
     'source/user/beamdyn/references.bib',
     'source/user/extptfm/bibliography.bib',
     'source/user/fast.farm/bibliography.bib',

docs/source/user/aerodyn/bibliography.bib

@@ -116,3 +116,10 @@ @article{ad-hammam2022
     volume= 1,
     number=1
 }
+
+@techreport{ad-hammam_NREL:2023,
+  title={Modeling the Yaw Behavior of Tail Fins for Small Wind Turbines: November 22, 2021-May 21, 2024},
+  author={Hammam, Mohamed M and Wood, David and Summerville, Brent},
+  year={2023},
+  institution={National Renewable Energy Laboratory (NREL), Golden, CO (United States)}
+}

docs/source/user/aerodyn/input.rst

@@ -883,22 +883,26 @@ An example of tail fin input file is given below:
     0         TFinIndMod  - Model for induced velocity calculation {0: none, 1:rotor-average} (switch)
     ====== Polar-based model ================================ [used only when TFinMod=1] 
     1        TFinAFID - Index of Tail fin airfoil number [1 to NumAFfiles]
-    ====== Unsteady slender body model  ===================== [used only when TFinMod=2] 
-    [TODO inputs for model 2]
+    ====== Unsteady slender body model  ===================== [used only when TFinMod=2]
+    0.9           TFinKp        - Tail fin moment of area about reference point
+    0.3,0.1,0.1   TFinSigma     - Tail fin empirical constant for vortex separation functions
+    40,60,60      TFinAStar     - Tail fin initial angles for vortex separation functions (deg)
+    3.1416        TFinKv        - Tail fin vortex lift coefficient
+    1.3           TFinCDc       - Tail fin drag coefficient
 
 General inputs
 ~~~~~~~~~~~~~~
 
-**TFinMod** Switch to select a model for the tail fin aerodynamics:
+``TFinMod`` is a switch to select a model for the tail fin aerodynamics:
 0) none (the aerodynamic forces are zero), 1) polar-based, 2) USB-based (see :numref:`TF-aerotheory`).
 (switch)
 
-**TFinArea** Area of the tail fin. (m^2)
+``TFinArea`` is the area of the tail fin. (m^2)
 This is the plan form area of the tail fin plate used to relate the local dynamic pressure and airfoil
 coefficients to aerodynamic loads. This value must not be negative and is only used when
 TFinMod is set to 1. (m^2)
 
-**TFinRefP_n** Undeflected position (:math:`x_{\text{ref},x_n},x_{\text{ref},y_n}, x_{\text{ref},z_n}`) of the tail fin from the tower top in nacelle coordinates.
+``TFinRefP_n`` is the undeflected position (:math:`x_{\text{ref},x_n},x_{\text{ref},y_n}, x_{\text{ref},z_n}`) of the tail fin from the tower top in nacelle coordinates.
 (formerly defined using ``TFinCPxn``,  ``TFinCPyn``, ``TFinCPzn``). 
 The distances defines the configuration for a furl angle of zero.
 For a typical upwind wind turbine, 
@@ -908,7 +912,7 @@ For a typical upwind wind turbine,
 See :numref:`figTFGeom` and :numref:`figTFcoord1`.
 (m)
 
-**TFinAngles** Angles (:math:`\theta_\text{skew},\theta_\text{tilt}, \theta_\text{bank}`) of the tail fin 
+``TFinAngles`` are the angles (:math:`\theta_\text{skew},\theta_\text{tilt}, \theta_\text{bank}`) of the tail fin 
 (formerly defined as ``TFinSkew``, ``TFinTilt``, ``TFinBank``).
 See :numref:`figTFGeom` and :numref:`figTFcoord1`. 
 These angles define the chordline at a furl angle of zero, where the chordline is assumed to be passing through the reference point.
@@ -925,7 +929,7 @@ This value must be greater than -180 and less than or equal to 180 degrees.
 
 
 
-**TFinIndMod**
+``TFinIndMod``
 Switch to select a model for the calculation of the velocity induced by the rotor and its wake on the tailfin (not the induced velocity from the tailfin wing).
 The options available are: 
 0) none (the induced velocity is zero)
@@ -936,7 +940,7 @@ The options available are:
 Polar-based model inputs
 ~~~~~~~~~~~~~~~~~~~~~~~~
 
-**TFinAFID**
+``TFinAFID``
 This integer tells AeroDyn which of the input airfoil files (``AFNames``) is assigned to the tail fin. For
 instance, a value of 2 means that the tail fin will use ``AFNames(2)`` for the local tail fin airfoil. 
 This value must be
@@ -945,7 +949,21 @@ between 1 and ``NumAFfiles`` and is only used when TFinMod is set to 1. (-)
 
 Unsteady slender body (USB) model inputs
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+Refer to :numref:`TF-aerotheory` and :cite:`ad-hammam_NREL:2023` for guidance on how to select parameters for the unsteady slender body theory based model.
 
-This option is currently not available and will be documented in a future release.
+``TFinKp``
+Potential lift coefficient for unsteady aerodynamics. ``TFinKp`` is used to calculate the potential flow contribution to the unsteady aerodynamic force on the tail fin.
 
+``TFinSigma``
+Tail fin empirical constants characterizing the decay of separation functions used in the unsteady aerodynamic model. The separation functions and their dependence on ``TFinSigma`` are described in :numref:`TF-aerotheory`.
 
+
+``TFinAStar``
+Tail fin characteristics angles for separation functions used in the unsteady aerodynamic model. The separation functions and their dependence on ``TFinAStar`` are described in :numref:`TF-aerotheory`.
+
+
+``TFinKv``
+Vortex lift coefficient for unsteady aerodynamics. ``TFinKv`` is used to calculate the vortex flow contribution to the unsteady aerodynamic force on the tail fin.
+
+``TFinCDc``
+Tail fin drag coefficient used for unsteady aerodynamic model. The drag on the tail fin significantly contributes to the normal force at high yaw angles.

docs/source/user/aerodyn/theory_tailfin.rst

@@ -150,8 +150,6 @@ Where :math:`\boldsymbol{V}_{\text{ind},\text{blade}}[i_b, i_r]` is the induced
 More advanced models could set the induced velocity to zero when outside of the wake boundary, or include a tower-shadow-like wake model. Such option is not yet available.
 
 
-
-
 Polar-based model
 -----------------
 
@@ -163,9 +161,25 @@ The user only needs to indicate the index `TFinAFIndex` within the list `AFNames
 Unsteady slender body model
 ---------------------------
 
-The unsteady slender body (USB) model is documented in :cite:`ad-hammam2022`.
+The unsteady aerodynamics of the tail fin is modeled based on Unsteady Slender Body Theory.
+The theory is extended to include the effect of high yaw angle :cite:`ad-hammam_NREL:2023`.
+
+The normal force on the tail fin can be described as the sum of three contributions (potential lift, vortex lift, and drag), weighted by separation functions :math:`x_i` as:
+
+.. math::  :label: TFUSBForce
+
+    N = \frac{\rho}{2} A_{tf} \bigg(  K_p x_1 V_{\text{rel},x} V_{\text{rel},y} +  \Big[x_2 K_v+(1- x_3)C_{Dc} \Big] V_{\text{rel},y}\big|V_{\text{rel},y}\big|\bigg)
+
+
+where :math:`\rho` is the density of air, :math:`A_{tf}` is the tail fin area, :math:`K_p` is the potential lift coefficient and :math:`K_v` is the vortex lift coefficient, and :math:`C_{Dc}` is the drag coefficient.
+:math:`x_i` are the separation functions calculated using a quasi-steady approximation as:
+
+.. math::  :label: TFUSBxiEquation
 
-The theory will be implemented and documented in a future release.
+    x_i = (1+exp{[\sigma_i (|\gamma_{tf}|-\alpha^*_i)]})^{-1}
 
 
+where :math:`\sigma_i` are empirical constants characterizing the decay of separation functions, :math:`\gamma_{tf}` is the yaw angle of the tail fin with respect to the free-stream wind (:math:`V_{\text{wind}}`), :math:`\alpha^*_i` are the characteristics angles for separation functions.
+:math:`x_i` takes on a value between 0 and 1, and are used to activate or deactivate a the contribution of potential lift, vortex lift and draft to the normal force on the tail fin.
 
+The normal force is assumed to act at the user defined reference point on the tail fin and the moment of the normal force is calculated accordingly.

docs/source/user/elastodyn/bibliography.bib

@@ -0,0 +1,6 @@
+@techreport{ed-hammam2023,
+  title={Modeling the Yaw Behavior of Tail Fins for Small Wind Turbines: November 22, 2021-May 21, 2024},
+  author={Hammam, Mohamed M and Wood, David and Summerville, Brent},
+  year={2023},
+  institution={National Renewable Energy Laboratory (NREL), Golden, CO (United States)}
+}

docs/source/user/elastodyn/index.rst

@@ -58,3 +58,4 @@ equations of FAST v7 and the ElastoDyn module of FAST v8 and OpenFAST.
    coordsys.rst
    input.rst
    theory.rst
+   zrefs.rst

docs/source/user/elastodyn/input.rst

@@ -231,6 +231,17 @@ Rotor-Teeter
 **TeetHSSp**    - Rotor-teeter hard-stop linear-spring constant (N-m/rad) [used only for 2 blades and when TeetMod=1]
 
 
+Yaw-Friction
+~~~~~~~~~~~~
+
+**YawFrctMod**  - Yaw-friction model {0: none, 1: friction without Fz term at the yaw bearing, 2: friction includes Fz term at yaw bearing, 3: user defined model} 
+
+**M_CSmax**     - Maximum Coulomb friction torque (N-m)[mu_s*D_eff when YawFrctMod=1 and Fz*mu_s*D_eff when YawFrctMod=2]
+
+**M_CD**        - Dynamic friction moment at null yaw rate (N-m) [mu_d*D_eff when YawFrctMod=1 and Fz*mu_d*D_eff when YawFrctMod=2]
+
+**sig_v**       - Viscous friction coefficient (N-m/(rad/s))
+
 
 Drivetrain
 ~~~~~~~~~~

docs/source/user/elastodyn/theory.rst

@@ -159,7 +159,41 @@ The total moment on the given degree of freedom is:
 
 
 
+.. _ed_yawfriction_theory:
 
+Yaw-friction
+------------
+A yaw-friction model is implemented in ElastoDyn based on a Coulomb-viscous approach.
+The yaw-friction moment as a function of yaw rate (:math:`\omega`) is shown below in :numref:`figYawFriction`
+
+.. _figYawFriction:
+.. figure:: figs/YawFrictionModel.jpg
+   :width: 60%
+
+   Yaw-friction model
+
+The yaw-friction torque :math:`M_f` can be calcualated as follows.
+
+if :math:`\omega\neq0`:
+
+.. math::
+   M_f = \mu_d\bar{D}\times\text{min}(0,F_z)\times\text{sign}(\omega) - \sigma_v\times\omega  
+ 
+if :math:`\omega=0` and :math:`\dot{\omega}\neq 0`:
+
+.. math::
+   M_f = \mu_d\bar{D}\times\text{min}(0,F_z)\times\text{sign}(\dot{\omega})
+
+if :math:`\omega=0` and :math:`\dot{\omega}=0`:
+
+.. math::
+   M_f = -\text{min}(\mu_s\bar{D}\times|\text{min}(0,F_z)|,|M_z|)\times\text{sign}(M_z)
+
+
+where :math:`\bar{D}` is the effective yaw-bearing race diameter, :math:`\mu_d` is the dynamic friction coefficient, :math:`\mu_s` is the static friction coefficient, :math:`F_z` is effective axial load on yaw-bearing, :math:`M-z` is the external torque on yaw-bearing.
+The static 'stiction' (where the static contribution exceeds the dynamic Coulomb friction) is only applied if both the yaw rotational velocity and acceleration at the current time-step are zero.
+The static portion of the friction is omitted if the rotational acceleration is not null, (sign(0) is taken as 1).
+This is to account for the fact that a 'warm' joint may not feel stiction when crossing through zero velocity in a dynamic sense :cite:`ed-hammam2023`
 
 
 .. _ed_dev_notes:

docs/source/user/elastodyn/zrefs.rst

@@ -0,0 +1,9 @@
+.. only:: html
+
+   References
+   ----------
+
+.. bibliography:: bibliography.bib
+   :labelprefix: ed-
+
+

modules/aerodyn/src/AeroDyn.f90

@@ -392,6 +392,11 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut
       p%rotors(iR)%TFin%TFinArea    = InputFileData%rotors(iR)%TFin%TFinArea
       p%rotors(iR)%TFin%TFinIndMod  = InputFileData%rotors(iR)%TFin%TFinIndMod
       p%rotors(iR)%TFin%TFinAFID    = InputFileData%rotors(iR)%TFin%TFinAFID
+      p%rotors(iR)%TFin%TFinKp      = InputFileData%rotors(iR)%TFin%TFinKp
+      p%rotors(iR)%TFin%TFinSigma   = InputFileData%rotors(iR)%TFin%TFinSigma
+      p%rotors(iR)%TFin%TFinAStar   = InputFileData%rotors(iR)%TFin%TFinAStar
+      p%rotors(iR)%TFin%TFinKv      = InputFileData%rotors(iR)%TFin%TFinKv
+      p%rotors(iR)%TFin%TFinCDc     = InputFileData%rotors(iR)%TFin%TFinCDc
    enddo
 
       !............................................................................................
@@ -4357,9 +4362,12 @@ SUBROUTINE TFin_CalcOutput(p, p_AD, u, m, y, ErrStat, ErrMsg )
    real(ReKi)              :: V_wnd(3)          ! wind velocity
    real(ReKi)              :: V_ind(3)          ! induced velocity
    real(ReKi)              :: V_str(3)          ! structural velocity
+   real(ReKi)              :: V_wnd_tf(3)          ! wind velocity
    real(ReKi)              :: force_tf(3)       ! force in tf system
    real(ReKi)              :: moment_tf(3)      ! moment in tf system
    real(ReKi)              :: alpha, Re, Cx, Cy, q ! Cl, Cd, Cm, 
+   real(ReKi)              :: x1, x2, x3,gamma_tf! scaling functions, gamma for unsteady modeling
+
    type(AFI_OutputType)    :: AFI_interp  ! Resulting values from lookup table
    integer(intKi)          :: ErrStat2
    character(ErrMsgLen)    :: ErrMsg2
@@ -4375,46 +4383,62 @@ SUBROUTINE TFin_CalcOutput(p, p_AD, u, m, y, ErrStat, ErrMsg )
 
    if (p%TFin%TFinIndMod==TFinIndMod_none) then
       V_ind = 0.0_ReKi
+
    elseif(p%TFin%TFinIndMod==TFinIndMod_rotavg) then
       ! TODO TODO
       print*,'TODO TailFin: compute rotor average induced velocity'
       V_ind = 0.0_ReKi 
+
    else
-      STOP ! Will never happen
+      call setErrStat(ErrID_Fatal, 'TailFin model unsupported', ErrStat, ErrMsg, 'TFin_CalcOutput')
+
    endif
-   V_rel       = V_wnd - V_str + V_ind
+
+   V_rel       = V_wnd - V_str + V_ind                          ! relative wind on tail fin
    V_rel_tf    = matmul(u%TFinMotion%Orientation(:,:,1), V_rel) ! from inertial to tf system
-   alpha       = atan2( V_rel_tf(2), V_rel_tf(1))               ! angle of attack
+   alpha       = atan2(V_rel_tf(2), V_rel_tf(1))                ! angle of attack
+   v_wnd_tf    = matmul(u%TFinMotion%Orientation(:,:,1), V_wnd) ! only used for calculation of x1,x2,x3
+   gamma_tf = atan2(v_wnd_tf(2), v_wnd_tf(1))                   ! only used for calculation of x1,x2,x3
    V_rel_orth2 = V_rel_tf(1)**2 + V_rel_tf(2)**2                ! square norm of Vrel in tf system
 
+   ! Initialize the tail fin forces to zero
+   force_tf(:)    = 0.0_ReKi
+   moment_tf(:)   = 0.0_ReKi
+
    if (p%TFin%TFinMod==TFinAero_none) then
-      y%TFinLoad%Force(1:3,1)  = 0.0_ReKi
-      y%TFinLoad%Moment(1:3,1) = 0.0_ReKi
+      ! Do nothing
 
    elseif (p%TFin%TFinMod==TFinAero_polar) then
-      ! Airfoil coefficients
+      ! Airfoil coefficients based model
       Re  = sqrt(V_rel_orth2) * p%TFin%TFinChord/p%KinVisc
       call AFI_ComputeAirfoilCoefs( alpha, Re, 0.0_ReKi,  p_AD%AFI(p%TFin%TFinAFID), AFI_interp, ErrStat2, ErrMsg2)
       call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
       Cx = -AFI_interp%Cl * sin(alpha) + AFI_interp%Cd * cos(alpha)
       Cy =  AFI_interp%Cl * cos(alpha) + AFI_interp%Cd * sin(alpha)
       ! Forces in tailfin system
       q = 0.5 * p%airDens * V_rel_orth2 * p%TFin%TFinArea
-      force_tf(:)    = 0.0_ReKi
-      moment_tf(:)    = 0.0_ReKi
+
       force_tf(1)    = Cx * q
       force_tf(2)    = Cy * q
-      force_tf(3)    = 0.0_ReKi
-      moment_tf(1:2) = 0.0_ReKi
       moment_tf(3)   = AFI_interp%Cm * q * p%TFin%TFinChord
-      ! Transfer to global
-      y%TFinLoad%Force(1:3,1)  = matmul(transpose(u%TFinMotion%Orientation(:,:,1)), force_tf)
-      y%TFinLoad%Moment(1:3,1) = matmul(transpose(u%TFinMotion%Orientation(:,:,1)), moment_tf)
 
    elseif (p%TFin%TFinMod==TFinAero_USB) then
-      call SetErrStat(ErrID_Fatal, 'Tail fin USB model not yet available', ErrStat, ErrMsg, RoutineName )
-      return
+      ! Unsteady aerodynamic model
+
+      ! Calculate separation function (quasi-steady)
+      x1 = 1.0_Reki/(1.0_Reki+exp(p%TFin%TFinSigma(1)*((ABS(gamma_tf)*180.0_ReKi/pi)-p%TFin%TFinAStar(1)))) 
+      x2 = 1.0_Reki/(1.0_Reki+exp(p%TFin%TFinSigma(2)*((ABS(gamma_tf)*180.0_ReKi/pi)-p%TFin%TFinAStar(2)))) 
+      x3 = 1.0_Reki/(1.0_Reki+exp(p%TFin%TFinSigma(3)*((ABS(gamma_tf)*180.0_ReKi/pi)-p%TFin%TFinAStar(3))))
+
+      ! Calculate unsteady force on tail fin
+      force_tf(2) = 0.5_ReKi * p%AirDens * p%TFin%TFinArea * &
+         (p%TFin%TFinKp * x1 * V_rel_tf(1) * V_rel_tf(2) + &
+         (x2 * p%TFin%TFinKv + (1-x3)*p%TFin%TFinCDc) * V_rel_tf(2) * ABS(V_rel_tf(2)))
    endif
+
+   ! Transfer to global
+   y%TFinLoad%Force(1:3,1)  = matmul(transpose(u%TFinMotion%Orientation(:,:,1)), force_tf)
+   y%TFinLoad%Moment(1:3,1) = matmul(transpose(u%TFinMotion%Orientation(:,:,1)), moment_tf)
 
    ! --- Store
    m%TFinAlpha  = alpha
@@ -4428,6 +4452,7 @@ SUBROUTINE TFin_CalcOutput(p, p_AD, u, m, y, ErrStat, ErrMsg )
    m%TFinM_i    = y%TFinLoad%Moment(1:3,1)
 
 END SUBROUTINE TFin_CalcOutput
+
 !----------------------------------------------------------------------------------------------------------------------------------
 !> This subroutine calculates the tower loads for the AeroDyn TowerLoad output mesh.
 SUBROUTINE ADTwr_CalcOutput(p, u, m, y, ErrStat, ErrMsg )

modules/aerodyn/src/AeroDyn_IO.f90

@@ -1266,15 +1266,20 @@ SUBROUTINE ReadTailFinInputs(FileName, TFData, UnEc, ErrStat, ErrMsg)
    call ParseVar(FileInfo_In, iLine, 'TFinAFID'  , TFData%TFinAFID      , ErrStat2, ErrMsg2, UnEc); if (Failed()) return;
    !====== Unsteady slender body model ===================== [used only when TFinMod=2]
    call ParseCom(FileInfo_in, iLine, DummyLine                          , ErrStat2, ErrMsg2, UnEc); if (Failed()) return;
+   call ParseVar(FileInfo_In, iLine, 'TFinKp'  , TFData%TFinKp      , ErrStat2, ErrMsg2, UnEc); if (Failed()) return;
+   call ParseAry(FileInfo_In, iLine, 'TFinSigma'  , TFData%TFinSigma, 3 , ErrStat2, ErrMsg2, UnEc); if (Failed()) return;
+   call ParseAry(FileInfo_In, iLine, 'TFinAStar', TFData%TFinAStar, 3 , ErrStat2, ErrMsg2, UnEc); if (Failed()) return;
+   call ParseVar(FileInfo_In, iLine, 'TFinKv'    , TFData%TFinKv       , ErrStat2, ErrMsg2, UnEc); if (Failed()) return;
+   call ParseVar(FileInfo_In, iLine, 'TFinCDc'   , TFData%TFinCDc       , ErrStat2, ErrMsg2, UnEc); if (Failed()) return;
+
    ! TODO
 
    ! --- Triggers
    TFData%TFinAngles = TFData%TFinAngles*D2R ! deg2rad
 
    ! --- Validation on the fly
-   !if (all((/TFinAero_none,TFinAero_polar, TFinAero_USB/) /= TFData%TFinMod)) then
-   if (all((/TFinAero_none,TFinAero_polar/) /= TFData%TFinMod)) then
-      call Fatal('TFinMod needs to be 0, or 1')
+   if (all((/TFinAero_none,TFinAero_polar,TFinAero_USB/) /= TFData%TFinMod)) then
+      call Fatal('TFinMod needs to be 0, 1 or 2')
    endif
    !if (all((/TFinIndMod_none,TFinIndMod_rotavg/) /= TFData%TFinIndMod)) then
    if (all((/TFinIndMod_none/) /= TFData%TFinIndMod)) then