Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add unsteady aerodynamic model for turbine tail fin #1874

Merged
merged 9 commits into from
Mar 11, 2024
7 changes: 7 additions & 0 deletions docs/source/user/aerodyn/bibliography.bib
Original file line number Diff line number Diff line change
Expand Up @@ -108,3 +108,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)}
}
36 changes: 27 additions & 9 deletions docs/source/user/aerodyn/input.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggest adding units after each variable.

For unitless variables, I would use (-). According to the theory/code, it looks like TFinSigma has the units of 1/degrees.

0.9 TFinKp - Tail fin moment of area about reference point
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This comment for TFinKp is not correct. TFinKp is the potential lift coefficient.

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,
Expand All @@ -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.
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wasn't there supposed to be preprocessor developed by Hammam et al so that these parameters could be set directly based on the shape of the tail fin (delta, rectangular, elliptic, etc.)? Was this ever developed? If so, could we share a link?


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.
22 changes: 18 additions & 4 deletions docs/source/user/aerodyn/theory_tailfin.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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
-----------------

Expand All @@ -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.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggest summarizing the limitations of the implementation here, i.e., tail fin chord assumed much smaller than tail fin arm length, and quasi-steady conditions (because your separation functions are implemented only quasi-statically). Perhaps there are others stated in Hammam et al.

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)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is a negative sign out front in Hammam et al. I presume that is not needed here due to the sign convention used in AeroDyn? I haven't confirmed this, but it may be good to state that.

Presumably this implementation matches results from Hammam et al?



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}
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't see an ABS() in Hammam et al. It makes sense to me that this would be needed, but I suggest mentioning that difference.



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.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why are the separation functions computed using ambient wind while the force calculation uses relative wind (including induction)? Shouldn't the same velocities (which I assume should be the relative velocity) be used for all?

: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.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Change "draft" to "drag here.


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.
55 changes: 40 additions & 15 deletions modules/aerodyn/src/AeroDyn.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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

!............................................................................................
Expand Down Expand Up @@ -4356,9 +4361,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
Expand All @@ -4374,46 +4382,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))))
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could use NWTC Library parameter R2D instead of 180/pi here.

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
Expand All @@ -4427,6 +4451,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 )
Expand Down
11 changes: 8 additions & 3 deletions modules/aerodyn/src/AeroDyn_IO.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
10 changes: 10 additions & 0 deletions modules/aerodyn/src/AeroDyn_Registry.txt
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,11 @@ typedef ^ TFinParameterType ReKi TFinChord - - - "Tail fin chord [u
typedef ^ TFinParameterType ReKi TFinArea - - - "Tail fin planform area [used only when TFinMod=1]" m^2
typedef ^ TFinParameterType IntKi TFinIndMod - - - "Model for induced velocity calculation {0=none, 1=rotor-average}" (switch)
typedef ^ TFinParameterType IntKi TFinAFID - - - "Index of Tail fin airfoil number [1 to NumAFfiles]" -
typedef ^ TFinParameterType ReKi TFinKp - - - "Tail fin potential lift coefficient for unsteady aerodynamics [used only when TFMod=2]" -
typedef ^ TFinParameterType ReKi TFinSigma 3 - - "Tail fin empirical constants characterizing the decay of separation functions [used only when TFMod=2]" -
typedef ^ TFinParameterType ReKi TFinAStar 3 - - "Tail fin characteristics angles for separation functions [used only when TFMod=2]" deg
typedef ^ TFinParameterType ReKi TFinKv - - - "Tail fin vortex lift coefficient for unsteady aerodynamics [used only when TFMod=2]" -
typedef ^ TFinParameterType ReKi TFinCDc - - - "Tail fin drag coefficient for unsteady aerodynamics [used only when TFMod=2]" -

# Tail Fin input file
typedef ^ TFinInputFileType IntKi TFinMod - - 0 "Tail fin aerodynamics model {0=none, 1=polar-based, 2=USB-based}" (switch)
Expand All @@ -63,6 +68,11 @@ typedef ^ TFinInputFileType ReKi TFinRefP_n 3 - - "Undeflected posit
typedef ^ TFinInputFileType ReKi TFinAngles 3 - - "Tail fin chordline skew, tilt, and bank angles about the reference point" (deg)
typedef ^ TFinInputFileType IntKi TFinIndMod - - - "Model for induced velocity calculation {0=none, 1=rotor-average}" (switch)
typedef ^ TFinInputFileType IntKi TFinAFID - - - "Index of Tail fin airfoil number [1 to NumAFfiles]" -
typedef ^ TFinInputFileType ReKi TFinKp - - - "Tail fin potential lift coefficient for unsteady aerodynamics [used only when TFMod=2]" -
typedef ^ TFinInputFileType ReKi TFinSigma 3 - - "Tail fin empirical constants characterizing the decay of separation functions [used only when TFMod=2]" -
typedef ^ TFinInputFileType ReKi TFinAStar 3 - - "Tail fin characteristics angles for separation functions [used only when TFMod=2]" deg
typedef ^ TFinInputFileType ReKi TFinKv - - - "Tail fin vortex lift coefficient for unsteady aerodynamics [used only when TFMod=2]" -
typedef ^ TFinInputFileType ReKi TFinCDc - - - "Tail fin drag coefficient for unsteady aerodynamics [used only when TFMod=2]" -



Expand Down
Loading
Loading