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

AD/AA: Add new TE definition, improve airfoil thickness calculation, simplify input #597

Merged
merged 21 commits into from
Feb 4, 2021
Merged
Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
21 commits
Select commit Hold shift + click to select a range
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
62 changes: 22 additions & 40 deletions docs/source/user/aerodyn-aeroacoustics/App-usage.rst
Original file line number Diff line number Diff line change
Expand Up @@ -97,12 +97,7 @@ models:
- **BluntMod** – Integer 0/1: flag to activate (**BluntMod=1**) the
trailing-edge bluntness – vortex shedding model, see :numref:`aa-TE-vortex`. If
the flag is set to 1, the trailing-edge geometry must be specified in
the file(s) listed in the field Blade Properties.

Next, the field Blade Properties lists three file names, often but not
necessarily identical, which contain the distributed properties
describing the detailed geometry of the trailing edge. These are
described in :numref:`aa-sec-TEgeom`.
the files as described in :numref:`aa-sec-BLinputs`.

The field Observer Locations contains the path to the file where the
number of observers (NrObsLoc) and the respective locations are
Expand Down Expand Up @@ -141,7 +136,7 @@ The file must be closed by an END command.

.. _aa-sec-BLinputs:

Boundary Layer Inputs
Boundary Layer Inputs and Trailing Edge Geometry
---------------------

When the flag **BLMod** is set equal to 2, pretabulated properties of the
Expand Down Expand Up @@ -188,6 +183,26 @@ outputs of XFoil. Because it is usually impossible to obtain these values for
the whole ranges of Reynolds numbers and angles of attack, the code is set to
adopt the last available values and print to screen a warning.

When the flag **BluntMod** is set to 1, the detailed geometry of the
trailing edge must also be defined along the span. Two inputs must be
provided, namely the angle, :math:`\Psi` between the suction and
pressure sides of the profile, right before the trailing-edge point, and
the height, :math:`h`, of the trailing edge. :math:`\Psi` must be
defined in degrees, while :math:`h` is in meters. Note that the BPM
trailing-edge bluntness model is very sensitive to these two parameters,
which, however, are often not easy to determine for real blades.
:numref:`aa-fig:GeomParamTE` shows the two inputs.

.. figure:: media/NoiseN011.png
:alt: Geometric parameters of the trailing-edge bluntness
:name: aa-fig:GeomParamTE
:width: 100.0%

Geometric parameters :math:`\mathbf{\Psi}` and
:math:`\mathbf{h}` of the trailing-edge bluntness

One value of :math:`\Psi` and one value of :math:`h` per file must be defined.
These values are not used if the flag **BluntMod** is set to 0.

.. container::
:name: aa-tab:AF20_BL
Expand Down Expand Up @@ -264,38 +279,5 @@ grid looks like the following:
:language: none


.. _aa-sec-TEgeom:

Trailing-Edge Geometry
----------------------

When the flag **BluntMod** is set to 1, the detailed geometry of the
trailing edge must be defined along the span. Two inputs must be
provided, namely the angle, :math:`\Psi,` between the suction and
pressure sides of the profile, right before the trailing-edge point, and
the height, :math:`h`, of the trailing edge. :math:`\Psi` must be
defined in degrees, while :math:`h` is in meters. Note that the BPM
trailing-edge bluntness model is very sensitive to these two parameters,
which, however, are often not easy to determine for real blades.
:numref:`aa-fig:GeomParamTE` shows the two inputs.

.. figure:: media/NoiseN011.png
:alt: Geometric parameters of the trailing-edge bluntness
:name: aa-fig:GeomParamTE
:width: 100.0%

Geometric parameters :math:`\mathbf{\Psi}` and
:math:`\mathbf{h}` of the trailing-edge bluntness

The two distributions must be defined with the same spanwise resolution
of the AeroDyn15 blade file, such as:

.. container::
:name: aa-tab:BladeProp

.. literalinclude:: example/BladeProp.dat
:linenos:
:language: none

.. [4]
https://github.com/OpenFAST/python-toolbox
138 changes: 71 additions & 67 deletions modules/aerodyn/src/AeroAcoustics.f90
Original file line number Diff line number Diff line change
Expand Up @@ -123,10 +123,10 @@ subroutine SetParameters( InitInp, InputFileData, p, ErrStat, ErrMsg )
CHARACTER(ErrMsgLen) :: ErrMsg2 ! temporary Error message if ErrStat / = ErrID_None
INTEGER(IntKi) :: ErrStat2 ! temporary Error status of the operation
INTEGER(IntKi) :: simcou,coun ! simple loop counter
INTEGER(IntKi) :: I,J,whichairfoil,K
INTEGER(IntKi) :: I,J,whichairfoil,K,i1_1,i10_1,i1_2,i10_2,iLE
character(*), parameter :: RoutineName = 'SetParameters'
LOGICAL :: tr,tri,exist
REAL(ReKi) :: val1,val2,f2,f4,lefttip,rightip,jumpreg
LOGICAL :: tr,tri,exist,LE_flag
REAL(ReKi) :: val1,val10,f2,f4,lefttip,rightip,jumpreg, dist1, dist10
! Initialize variables for this routine
ErrStat = ErrID_None
ErrMsg = ""
Expand Down Expand Up @@ -241,18 +241,19 @@ subroutine SetParameters( InitInp, InputFileData, p, ErrStat, ErrMsg )
call AllocAry(p%StallStart,p%NumBlNds,p%NumBlades,'p%StallStart',ErrStat2,ErrMsg2); if(Failed()) return
p%StallStart(:,:) = 0.0_ReKi

do i=1,p%NumBlades
p%TEThick(:,i) = InputFileData%BladeProps(i)%TEThick(:) !
p%TEAngle(:,i) = InputFileData%BladeProps(i)%TEAngle(:) !
do i=1,p%NumBlades
do j=1,p%NumBlNds
whichairfoil = p%BlAFID(j,i)
if(p%AFInfo(whichairfoil)%NumTabs /=1 ) then
call SetErrStat(ErrID_Fatal, 'Number of airfoil tables within airfoil file different than 1, which is not supported.', ErrStat2, ErrMsg2, RoutineName )
if(Failed()) return
endif
p%StallStart(j,i) = p%AFInfo(whichairfoil)%Table(1)%UA_BL%alpha1*180/PI ! approximate stall angle of attack [deg] (alpha1 in [rad])
whichairfoil = p%BlAFID(j,i)
p%TEThick(j,i) = InputFileData%BladeProps(whichairfoil)%TEThick
p%TEAngle(j,i) = InputFileData%BladeProps(whichairfoil)%TEAngle

if(p%AFInfo(whichairfoil)%NumTabs /=1 ) then
call SetErrStat(ErrID_Fatal, 'Number of airfoil tables within airfoil file different than 1, which is not supported.', ErrStat2, ErrMsg2, RoutineName )
if(Failed()) return
endif
p%StallStart(j,i) = p%AFInfo(whichairfoil)%Table(1)%UA_BL%alpha1*180/PI ! approximate stall angle of attack [deg] (alpha1 in [rad])
enddo
end do
enddo

call AllocAry(p%BlSpn, p%NumBlNds, p%NumBlades, 'p%BlSpn' , ErrStat2, ErrMsg2); if(Failed()) return
call AllocAry(p%BlChord, p%NumBlNds, p%NumBlades, 'p%BlChord', ErrStat2, ErrMsg2); if(Failed()) return
Expand Down Expand Up @@ -340,60 +341,63 @@ subroutine SetParameters( InitInp, InputFileData, p, ErrStat, ErrMsg )
if(Failed()) return
endif

! If simplified guidati is on, calculate the airfoil thickness from input airfoil coordinates
! If simplified guidati is on, calculate the airfoil thickness at 1% and at 10% chord from input airfoil coordinates
IF (p%IInflow .EQ. 2) THEN
! Calculate the Thickness @ 1% chord and @ 10% chord (normalized thickness)
call AllocAry(p%AFThickGuida,2,size(p%AFInfo), 'p%AFThickGuida', errStat2, errMsg2); if(Failed()) return
p%AFThickGuida=0.0_Reki

DO k=1,size(p%AFInfo) ! for each airfoil interpolation
tri=.true.;tr=.true.;
do i=2,size(p%AFInfo(k)%X_Coord)
if ( (p%AFInfo(k)%X_Coord(i)+p%AFInfo(k)%Y_Coord(i)) .eq. 0) then
!print*,i
goto 174
endif
if ( p%AFInfo(k)%X_Coord(i) .eq. 0.1) then
val1=p%AFInfo(k)%Y_Coord(i)
elseif ( (p%AFInfo(k)%X_Coord(i) .lt. 0.1) .and. (tri) ) then
val1=( abs(p%AFInfo(k)%X_Coord(i-1)-0.1)*p%AFInfo(k)%Y_Coord(i) + &
abs(p%AFInfo(k)%X_Coord(i)-0.1)*p%AFInfo(k)%Y_Coord(i-1))/ &
(abs(p%AFInfo(k)%X_Coord(i-1)-0.1)+abs(p%AFInfo(k)%X_Coord(i)-0.1))

tri=.false.
elseif (p%AFInfo(k)%X_Coord(i) .eq. 0.01) then
val2=p%AFInfo(k)%Y_Coord(i)
elseif ( (p%AFInfo(k)%X_Coord(i) .lt. 0.01) .and. (tr) ) then
val2=( abs(p%AFInfo(k)%X_Coord(i-1)-0.01)*p%AFInfo(k)%Y_Coord(i) + &
abs(p%AFInfo(k)%X_Coord(i)-0.01)*p%AFInfo(k)%Y_Coord(i-1))/ &
(abs(p%AFInfo(k)%X_Coord(i-1)-0.01)+abs(p%AFInfo(k)%X_Coord(i)-0.01))
tr=.false.
endif
enddo

174 tri=.true.;tr=.true.;
do j=i,size(p%AFInfo(k)%X_Coord)
if ( p%AFInfo(k)%X_Coord(j) .eq. 0.1) then
val1=abs(p%AFInfo(k)%Y_Coord(j)) + abs(val1)
elseif ( (p%AFInfo(k)%X_Coord(j) .gt. 0.1) .and. (tri) ) then
val1=abs(val1)+abs((abs(p%AFInfo(k)%X_Coord(j-1)-0.1)*p%AFInfo(k)%Y_Coord(j)+ &
abs(p%AFInfo(k)%X_Coord(j)-0.1)*p%AFInfo(k)%Y_Coord(j-1))/&
(abs(p%AFInfo(k)%X_Coord(j-1)-0.1)+abs(p%AFInfo(k)%X_Coord(j)-0.1)));
tri=.false.
elseif (p%AFInfo(k)%X_Coord(j) .eq. 0.01) then
val2=abs(p%AFInfo(k)%Y_Coord(j)) + abs(val2)
elseif ( (p%AFInfo(k)%X_Coord(j) .gt. 0.01) .and. (tr) ) then
val2=abs(val2)+abs((abs(p%AFInfo(k)%X_Coord(j-1)-0.01)*p%AFInfo(k)%Y_Coord(j)+ &
abs(p%AFInfo(k)%X_Coord(j)-0.01)*p%AFInfo(k)%Y_Coord(j-1))/&
(abs(p%AFInfo(k)%X_Coord(j-1)-0.01)+abs(p%AFInfo(k)%X_Coord(j)-0.01)));
tr=.false.
endif
enddo

p%AFThickGuida(1,k)=val2 ! 1 % chord thickness
p%AFThickGuida(2,k)=val1 ! 10 % chord thickness

! IF ((MIN(p%AFInfo(k)%X_Coord) < 0.) .or. (MAX(p%AFInfo(k)%X_Coord) > 0.)) THEN
! call SetErrStat ( ErrID_Fatal,'The coordinates of airfoil '//trim(num2lstr(k))//' are mot defined between x=0 and x=1. Code stops.' ,ErrStat, ErrMsg, RoutineName )
! ENDIF

! Flip the flag when LE is found and find index
LE_flag = .False.
DO i=3,size(p%AFInfo(k)%X_Coord)
IF (LE_flag .eqv. .False.) THEN
IF (p%AFInfo(k)%X_Coord(i) - p%AFInfo(k)%X_Coord(i-1) > 0.) THEN
LE_flag = .TRUE.
iLE = i
ENDIF
ENDIF
ENDDO

! From LE toward TE
dist1 = ABS( p%AFInfo(k)%X_Coord(iLE) - 0.01)
dist10 = ABS( p%AFInfo(k)%X_Coord(iLE) - 0.10)
DO i=iLE+1,size(p%AFInfo(k)%X_Coord)
IF (ABS(p%AFInfo(k)%X_Coord(i) - 0.01) < dist1) THEN
i1_1 = i
dist1 = ABS(p%AFInfo(k)%X_Coord(i) - 0.01)
ENDIF
IF (ABS(p%AFInfo(k)%X_Coord(i) - 0.1) < dist10) THEN
i10_1 = i
dist10 = ABS(p%AFInfo(k)%X_Coord(i) - 0.1)
ENDIF
ENDDO

! From TE to LE
dist1 = 0.99
dist10 = 0.90
DO i=1,iLE-1
IF (ABS(p%AFInfo(k)%X_Coord(i) - 0.01) < dist1) THEN
i1_2 = i
dist1 = ABS(p%AFInfo(k)%X_Coord(i) - 0.01)
ENDIF
IF (ABS(p%AFInfo(k)%X_Coord(i) - 0.1) < dist10) THEN
i10_2 = i
dist10 = ABS(p%AFInfo(k)%X_Coord(i) - 0.1)
ENDIF
ENDDO

val1 = p%AFInfo(k)%Y_Coord(i1_1) - p%AFInfo(k)%Y_Coord(i1_2)
val10 = p%AFInfo(k)%Y_Coord(i10_1) - p%AFInfo(k)%Y_Coord(i10_2)

p%AFThickGuida(1,k)=val1 ! 1 % chord thickness
p%AFThickGuida(2,k)=val10 ! 10 % chord thickness
ENDDO
ENDIF ! If simplified guidati is on, calculate the airfoil thickness
ENDIF

!! for turbulence intensity calculations on the fly every 5 meter the whole rotor area is divided vertically to store flow fields in each region
jumpreg=7
Expand Down Expand Up @@ -1063,8 +1067,8 @@ SUBROUTINE CalcAeroAcousticsOutput(u,p,m,xd,y,errStat,errMsg)
!--------Turbulent Boundary Layer Trailing Edge Noise----------------------------!
IF ( (p%ITURB .EQ. 1) .or. (p%ITURB .EQ. 2) ) THEN
CALL TBLTE(AlphaNoise,p%BlChord(J,I),UNoise,m%ChordAngleTE(K,J,I),m%SpanAngleTE(K,J,I), &
elementspan,m%rTEtoObserve(K,J,I), p, j,i,k,m%d99Var(2),m%dstarVar(1),m%dstarVar(2),p%StallStart(J,I), &
m%SPLP,m%SPLS,m%SPLALPH,m%SPLTBL,errStat2,errMsg2 )
elementspan,m%rTEtoObserve(K,J,I), p, j,i,k,m%d99Var(2),m%dstarVar(1),m%dstarVar(2),p%StallStart(J,I), &
m%SPLP,m%SPLS,m%SPLALPH,m%SPLTBL,errStat2,errMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
IF (p%ITURB .EQ. 2) THEN
m%SPLP=0.0_ReKi;m%SPLS=0.0_ReKi;m%SPLTBL=0.0_ReKi;
Expand All @@ -1078,9 +1082,9 @@ SUBROUTINE CalcAeroAcousticsOutput(u,p,m,xd,y,errStat,errMsg)
!--------Blunt Trailing Edge Noise----------------------------------------------!
IF ( p%IBLUNT .EQ. 1 ) THEN
CALL BLUNT(AlphaNoise,p%BlChord(J,I),UNoise,m%ChordAngleTE(K,J,I),m%SpanAngleTE(K,J,I), &
elementspan,m%rTEtoObserve(K,J,I),p%TEThick(J,I),p%TEAngle(J,I), &
p, m%d99Var(2),m%dstarVar(1),m%dstarVar(2),m%SPLBLUNT,p%StallStart(J,I),errStat2,errMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
elementspan,m%rTEtoObserve(K,J,I),p%TEThick(J,I),p%TEAngle(J,I), &
p, m%d99Var(2),m%dstarVar(1),m%dstarVar(2),m%SPLBLUNT,p%StallStart(J,I),errStat2,errMsg2 )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
ENDIF
!--------Tip Noise--------------------------------------------------------------!
IF ( (p%ITIP .EQ. 1) .AND. (J .EQ. p%NumBlNds) ) THEN
Expand Down Expand Up @@ -1445,7 +1449,7 @@ SUBROUTINE TBLTE(ALPSTAR,C,U,THETA,PHI,L,R,p,jj,ii,kk,d99Var2,dstarVar1,dstarVar

LOGICAL :: SWITCH !!LOGICAL FOR COMPUTATION OF ANGLE OF ATTACK CONTRIBUTION



ErrStat = ErrID_None
ErrMsg = ""
Expand Down
Loading