Skip to content

Commit

Permalink
Add unsteady aerodynamic model for turbine tail fin
Browse files Browse the repository at this point in the history
  • Loading branch information
@abhineet-gupta
abhineet-gupta committed Nov 15, 2023
1 parent ae56ab1 commit b4e3f4f
Show file tree
Hide file tree
Showing 6 changed files with 191 additions and 20 deletions.
7 changes: 7 additions & 0 deletions docs/source/user/aerodyn/bibliography.bib
View file
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)}
}
19 changes: 14 additions & 5 deletions docs/source/user/aerodyn/theory_tailfin.rst
View file
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 @@ -160,12 +158,23 @@ The tabulated data are provided as part of the list of airfoils given with `AFNa
The user only needs to indicate the index `TFinAFIndex` within the list `AFNames` to indicate which polar to use for the tail fin.


Unsteady slender body model
Unsteady 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

.. math:: :label: TFUSBForce

N = \frac{\rho}{2} A_{tf} K_p x_1 V_x V_y + \frac{\rho}{2} A_{tf} \Big[x_2 K_v+(1- x_3)C_{Dc} \Big] V_y|V_y|.


The theory will be implemented and documented in a future release.
And the moment on the tail fin about the apex can be described as:

.. math:: :label: TFUSBMoment

M_a = \frac{\rho}{2}A_{tf}x_{cp}x_1 K_p V_x V_y + \frac{\rho}{2}A_{tf}x_{cp}\Big[x_2K_v + (1-x_3)C_{Dc}\Big]V_y|V_y|

where :math:`A_{tf}` is the tail fin area, :math:`K_p` is the potential flow constant and :math:`K_v` is the vortex flow cosntant, :math:`x_i` are the separation function, and :math:`C_{Dc}` is the drag coefficient.
53 changes: 41 additions & 12 deletions modules/aerodyn/src/AeroDyn.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -392,6 +392,12 @@ 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%TFinCp = InputFileData%rotors(iR)%TFin%TFinCp
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 @@ -4358,7 +4364,10 @@ SUBROUTINE TFin_CalcOutput(p, p_AD, u, m, y, ErrStat, ErrMsg )
real(ReKi) :: V_str(3) ! structural 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) :: alpha, Re, Cx, Cy, q, tfingamma ! Cl, Cd, Cm,
! USB variables
real(ReKi) :: x1, x2, x3 ! scaling functions for different contributions on CN

type(AFI_OutputType) :: AFI_interp ! Resulting values from lookup table
integer(intKi) :: ErrStat2
character(ErrMsgLen) :: ErrMsg2
Expand All @@ -4382,13 +4391,18 @@ SUBROUTINE TFin_CalcOutput(p, p_AD, u, m, y, ErrStat, ErrMsg )
STOP ! Will never happen
endif
V_rel = V_wnd - V_str + V_ind
print *,'V_wnd'
print *,V_wnd
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
V_rel_orth2 = V_rel_tf(1)**2 + V_rel_tf(2)**2 ! square norm of Vrel in tf system


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
Expand All @@ -4399,21 +4413,35 @@ SUBROUTINE TFin_CalcOutput(p, p_AD, u, m, y, ErrStat, ErrMsg )
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
!Calculate separation functions
!x1 = 1.0_Reki/(1.0_Reki+exp(p%TFin%TFinSigma(1)*((ABS(alpha)*180.0_ReKi/pi)-p%TFin%TFinAStar(1))))
!x2 = 1.0_Reki/(1.0_Reki+exp(p%TFin%TFinSigma(2)*((ABS(alpha)*180.0_ReKi/pi)-p%TFin%TFinAStar(2))))
!x3 = 1.0_Reki/(1.0_Reki+exp(p%TFin%TFinSigma(3)*((ABS(alpha)*180.0_ReKi/pi)-p%TFin%TFinAStar(3))))

tfingamma = atan2(u%TFinMotion%orientation(2,1,1),u%TFinMotion%orientation(1,1,1))
x1 = 1.0_Reki/(1.0_Reki+exp(p%TFin%TFinSigma(1)*((ABS(tfingamma)*180.0_ReKi/pi)-p%TFin%TFinAStar(1))))
x2 = 1.0_Reki/(1.0_Reki+exp(p%TFin%TFinSigma(2)*((ABS(tfingamma)*180.0_ReKi/pi)-p%TFin%TFinAStar(2))))
x3 = 1.0_Reki/(1.0_Reki+exp(p%TFin%TFinSigma(3)*((ABS(tfingamma)*180.0_ReKi/pi)-p%TFin%TFinAStar(3))))

! print *,alpha*180.0_ReKi/pi
! print *,alpha

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)))
! moment_tf(3) = force_tf(2) * p%Tfin%TFinCp * p%TFin%TFinChord

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 +4455,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
12 changes: 9 additions & 3 deletions modules/aerodyn/src/AeroDyn_IO.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1266,15 +1266,21 @@ 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 ParseVar(FileInfo_In, iLine, 'TFinCp' , TFData%TFinCp , 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
12 changes: 12 additions & 0 deletions modules/aerodyn/src/AeroDyn_Registry.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,12 @@ 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 emperical constant for vortex separation functions [used only when TFMod=2]" -
typedef ^ TFinParameterType ReKi TFinCp - - - "Tail fin emperical constant for vortex separation functions [used only when TFMod=2]" -
typedef ^ TFinParameterType ReKi TFinSigma 3 - - "Tail fin emperical constants for vortex separation functions [used only when TFMod=2]" -
typedef ^ TFinParameterType ReKi TFinAStar 3 - - "Tail fin initial angles for vortex separation functions [used only when TFMod=2]" deg
typedef ^ TFinParameterType ReKi TFinKv - - - "Tail fin vortex lift coefficient [used only when TFMod=2]" -
typedef ^ TFinParameterType ReKi TFinCDc - - - "Tail fin drag coefficient [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 +69,12 @@ 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 emperical constant for vortex separation functions [used only when TFMod=2]" -
typedef ^ TFinInputFileType ReKi TFinCp - - - "Tail fin emperical constant for vortex separation functions [used only when TFMod=2]" -
typedef ^ TFinInputFileType ReKi TFinSigma 3 - - "Tail fin emperical constants for vortex separation functions [used only when TFMod=2]" -
typedef ^ TFinInputFileType ReKi TFinAStar 3 - - "Tail fin initial angles for vortex separation functions [used only when TFMod=2]" deg
typedef ^ TFinInputFileType ReKi TFinKv - - - "Tail fin vortex lift coefficient [used only when TFMod=2]" -
typedef ^ TFinInputFileType ReKi TFinCDc - - - "Tail fin drag coefficient [used only when TFMod=2]" -



Expand Down
Loading

0 comments on commit b4e3f4f

Please sign in to comment.