Skip to content

Commit

Permalink
OLAF: Bug Fix twist not computed and wrong Fn (OpenFAST#2348)
Browse files Browse the repository at this point in the history
  • Loading branch information
@ebranlard
ebranlard committed Jul 27, 2024
1 parent 4b8eba6 commit 4430e82
Show file tree
Hide file tree
Showing 5 changed files with 93 additions and 62 deletions.
99 changes: 56 additions & 43 deletions modules/aerodyn/src/AeroDyn.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -694,7 +694,9 @@ subroutine Init_MiscVars(m, p, p_AD, u, y, errStat, errMsg)
call AllocAry( m%TwrClrnc, p%NumBlNds, p%NumBlades, 'm%TwrClrnc', ErrStat2, ErrMsg2 )
call SetErrStat( errStat2, errMsg2, errStat, errMsg, RoutineName )
end if
call AllocAry( m%Curve, p%NumBlNds, p%NumBlades, 'm%Curve', ErrStat2, ErrMsg2 )
call AllocAry( m%Cant, p%NumBlNds, p%NumBlades, 'm%Cant', ErrStat2, ErrMsg2 )
call SetErrStat( errStat2, errMsg2, errStat, errMsg, RoutineName )
call AllocAry( m%Toe, p%NumBlNds, p%NumBlades, 'm%Toe', ErrStat2, ErrMsg2 )
call SetErrStat( errStat2, errMsg2, errStat, errMsg, RoutineName )
call AllocAry( m%X, p%NumBlNds, p%NumBlades, 'm%X', ErrStat2, ErrMsg2 )
call SetErrStat( errStat2, errMsg2, errStat, errMsg, RoutineName )
Expand Down Expand Up @@ -3120,18 +3122,20 @@ subroutine SetInputsForBEMT(p, p_AD, u, RotInflow, m, indx, errStat, errMsg)
if (ErrStat >= AbortErrLev) return

!..........................
! Set main geometry parameters (orientatioAnnulus, Curve, rLocal)
! Set main geometry parameters (orientatioAnnulus, Twist, Toe, Cant, rLocal)
!..........................
! TODO (EB): For harmonization between BEM and OLAF we should always compute R_li, r_Local, Twist, Toe, Cant
! BEM would then switch below between an "orientationMomentum", either Annulus (R_li) or NoPitchSweepPitch (R_wi)
if (p%AeroProjMod==APM_BEM_NoSweepPitchTwist .or. p%AeroProjMod==APM_LiftingLine) then

! orientationAnnulus and curve
if (p%AeroProjMod==APM_BEM_NoSweepPitchTwist) then
call Calculate_MeshOrientation_NoSweepPitchTwist(p, u, m, ErrStat=ErrStat, ErrMsg=ErrMsg, thetaBladeNds=thetaBladeNds)
call Calculate_MeshOrientation_NoSweepPitchTwist(p, u, m, thetaBladeNds, m%Toe, m%Cant, ErrStat=ErrStat, ErrMsg=ErrMsg)
else
call Calculate_MeshOrientation_LiftingLine(p, u, m, ErrStat=ErrStat, ErrMsg=ErrMsg, thetaBladeNds=thetaBladeNds)
call Calculate_MeshOrientation_LiftingLine(p, u, m, thetaBladeNds, m%Toe, m%Cant, ErrStat=ErrStat, ErrMsg=ErrMsg)
endif

! local radius (normalized distance from rotor centerline)
! local radius (normalized distance from rotor centerline) NOTE: unfortunate calculation, see comment above for harmonization
do k=1,p%NumBlades
call Calculate_MeshOrientation_Rel2Hub(u%BladeMotion(k), u%HubMotion, x_hat_disk, elemPosRelToHub_save=elemPosRelToHub, elemPosRotorProj_save=elemPosRotorProj)
do j=1,p%NumBlNds
Expand All @@ -3144,6 +3148,8 @@ subroutine SetInputsForBEMT(p, p_AD, u, RotInflow, m, indx, errStat, errMsg)

! Determine current azimuth angle and pitch axis vector of blade k, element j
call Calculate_MeshOrientation_Rel2Hub(u%BladeMotion(k), u%HubMotion, x_hat_disk, m%orientationAnnulus(:,:,:,k), elemPosRelToHub_save=elemPosRelToHub, elemPosRotorProj_save=elemPosRotorProj)
! Twist (aero+elastic), Toe, Cant (instantaneous and local), include elastic deformation
call TwistToeCant_FromLocalPolar(u%BladeMotion(k), m%orientationAnnulus(:,:,:,k), thetaBladeNds, m%Toe(:,k), m%Cant(:,k))

!..........................
! Compute local radius
Expand Down Expand Up @@ -3177,10 +3183,10 @@ subroutine SetInputsForBEMT(p, p_AD, u, RotInflow, m, indx, errStat, errMsg)


!..........................
! local blade angles
! local blade angles passed to BEM
!..........................
if (p%AeroProjMod==APM_BEM_NoSweepPitchTwist .or. p%AeroProjMod==APM_LiftingLine) then
! Theta
! Local and instantaneous blade twist+pitch (aerodynamic + elastic), cant and toe (include elastic deformation)
do k=1,p%NumBlades
do j=1,p%NumBlNds
m%BEMT_u(indx)%theta(j,k) = thetaBladeNds(j,k) ! local pitch + twist (aerodyanmic + elastic) angle of the jth node in the kth blade
Expand All @@ -3193,16 +3199,9 @@ subroutine SetInputsForBEMT(p, p_AD, u, RotInflow, m, indx, errStat, errMsg)
elseif (p%AeroProjMod==APM_BEM_Polar) then
do k=1,p%NumBlades
do j=1,p%NumBlNds
! Get local blade cant angle and twist
orientation = matmul( u%BladeMotion(k)%Orientation(:,:,j), transpose( m%orientationAnnulus(:,:,j,k) ) )
theta = EulerExtract( orientation )
! Get toe angle
m%BEMT_u(indx)%toeAngle(j,k) = theta(1)
! cant angle (including aeroelastic deformation)
m%BEMT_u(indx)%cantAngle(j,k) = theta(2)
m%Curve(j,k) = theta(2)
! twist (including pitch and aeroelastic deformation)
m%BEMT_u(indx)%theta(j,k) = -theta(3)
m%BEMT_u(indx)%theta(j,k) = thetaBladeNds(j,k)
m%BEMT_u(indx)%toeAngle(j,k) = m%Toe(j,k)
m%BEMT_u(indx)%cantAngle(j,k) = m%Cant(j,k)
end do !j=nodes
end do !k=blades
else
Expand Down Expand Up @@ -3378,6 +3377,27 @@ subroutine StorePitchAndAzimuth(p, u, m, ErrStat,ErrMsg)
enddo

endsubroutine StorePitchAndAzimuth
!----------------------------------------------------------------------------------------------------------------------------------
!> Instantaneous and local Twist Toe Cant angles from local polar to section
!! Note: could also be placed in Calculate_MeshOrientation_Rel2Hub
subroutine TwistToeCant_FromLocalPolar(secMesh, R_li, twist, toe, cant)
type(MeshType), intent(in ) :: secMesh !< Blade section mesh "BladeMotion"
real(R8Ki), intent(in ) :: R_li(3,3,secMesh%NNodes) !< Orientation from inertial (i) to local polar (l), aka "orientationAnnulus"
real(R8Ki), intent(out ) :: twist(secMesh%NNodes) !< Twist
real(ReKi), intent(out ) :: toe (secMesh%NNodes) !< Toe
real(ReKi), intent(out ) :: cant (secMesh%NNodes) !< Cant
real(R8Ki) :: R_sl(3,3) !< Orientation from local polar to section
integer(intKi) :: j !< loop counter for nodes
real(R8Ki) :: thetas(3) !< Euler angles
do j = 1, secMesh%NNodes
R_sl = matmul( secMesh%Orientation(:,:,j), transpose( R_li(:,:,j) ) ) ! From local polar to section - R_sec_i R_i_annulus
thetas = EulerExtract( R_sl )
toe(j) = real( thetas(1), ReKi) ! toe angle
cant(j) = real( thetas(2), ReKi) ! cant angle (including aeroelastic deformation)
twist(j) = -thetas(3) ! twist (including pitch and aeroelastic deformation)
end do
end subroutine TwistToeCant_FromLocalPolar

!----------------------------------------------------------------------------------------------------------------------------------
subroutine Calculate_MeshOrientation_Rel2Hub(Mesh1, HubMotion, x_hat_disk, orientationAnnulus, elemPosRelToHub_save, elemPosRotorProj_save)
TYPE(MeshType), intent(in) :: Mesh1 !< either BladeMotion or BladeRootMotion mesh
Expand Down Expand Up @@ -3427,25 +3447,16 @@ subroutine Calculate_MeshOrientation_Rel2Hub(Mesh1, HubMotion, x_hat_disk, orien
if (present(elemPosRotorProj_save)) elemPosRotorProj_save(:,j) = elemPosRotorProj
end do

! orientation = matmul( Mesh1(k)%Orientation(:,:,j), transpose( orientationAnnulus(:,:,j) ) )
! theta = EulerExtract( orientation )
! ! Get toe angle
! toeAngle(j) = theta(1)
! ! cant angle (including aeroelastic deformation)
! cantAngle(j) = theta(2)
! Curve(j) = theta(2)
! ! twist (including pitch and aeroelastic deformation)
! thetaNds(j) = -theta(3)

end subroutine Calculate_MeshOrientation_Rel2Hub
!----------------------------------------------------------------------------------------------------------------------------------
! Calculate_MeshOrientation_NoSweepPitchTwist sets orientationAnnulus, Curve and potential Blades nodes angles
subroutine Calculate_MeshOrientation_NoSweepPitchTwist(p, u, m, thetaBladeNds, toeBladeNds, ErrStat, ErrMsg)
subroutine Calculate_MeshOrientation_NoSweepPitchTwist(p, u, m, twist, toe, cant, ErrStat, ErrMsg)
type(RotParameterType), intent(in ) :: p !< AD parameters
type(RotInputType), intent(in ) :: u !< AD Inputs at Time
type(RotMiscVarType), intent(inout) :: m !< Misc/optimization variables
real(R8Ki), optional, intent( out) :: thetaBladeNds(p%NumBlNds,p%NumBlades)
real(R8Ki), optional, intent( out) :: toeBladeNds(p%NumBlNds,p%NumBlades)
real(R8Ki), optional, intent( out) :: twist(p%NumBlNds,p%NumBlades)
real(ReKi), optional, intent( out) :: toe(p%NumBlNds,p%NumBlades)
real(ReKi), optional, intent( out) :: cant(p%NumBlNds,p%NumBlades)
integer(IntKi), intent( out) :: ErrStat !< Error status of the operation
character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None
real(R8Ki) :: theta(3)
Expand Down Expand Up @@ -3483,9 +3494,9 @@ subroutine Calculate_MeshOrientation_NoSweepPitchTwist(p, u, m, thetaBladeNds, t
call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
theta = EulerExtract( orientation ) !root(k)WithoutPitch_theta(j)_blade(k)

m%Curve( j,k) = theta(2) ! save value for possible output later
if (present(thetaBladeNds)) thetaBladeNds(j,k) = -theta(3) ! local pitch + twist (aerodyanmic + elastic) angle of the jth node in the kth blade
if (present(toeBladeNds )) toeBladeNds( j,k) = theta(1)
if (present(cant)) cant (j,k) = theta(2) ! save value for possible output later
if (present(twist)) twist(j,k) = -theta(3) ! local pitch + twist (aerodyanmic + elastic) angle of the jth node in the kth blade
if (present(toe )) toe( j,k) = theta(1)

theta(1) = 0.0_ReKi
theta(3) = 0.0_ReKi
Expand All @@ -3495,15 +3506,16 @@ subroutine Calculate_MeshOrientation_NoSweepPitchTwist(p, u, m, thetaBladeNds, t
end do !k=blades
end subroutine Calculate_MeshOrientation_NoSweepPitchTwist
!----------------------------------------------------------------------------------------------------------------------------------
subroutine Calculate_MeshOrientation_LiftingLine(p, u, m, thetaBladeNds, toeBladeNds, ErrStat, ErrMsg)
subroutine Calculate_MeshOrientation_LiftingLine(p, u, m, twist, toe, cant, ErrStat, ErrMsg)
type(RotParameterType), intent(in ) :: p !< AD parameters
type(RotInputType), intent(in ) :: u !< AD Inputs at Time
type(RotMiscVarType), intent(inout) :: m !< Misc/optimization variables
real(R8Ki), optional, intent( out) :: thetaBladeNds(p%NumBlNds,p%NumBlades)
real(R8Ki), optional, intent( out) :: toeBladeNds(p%NumBlNds,p%NumBlades)
real(R8Ki), intent( out) :: twist(p%NumBlNds,p%NumBlades)
real(ReKi), intent( out) :: toe(p%NumBlNds,p%NumBlades)
real(ReKi), intent( out) :: cant(p%NumBlNds,p%NumBlades)
integer(IntKi), intent( out) :: ErrStat !< Error status of the operation
character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None
real(R8Ki) :: theta(3)
real(R8Ki) :: thetas(3)
real(R8Ki) :: orientation(3,3)
integer(intKi) :: j ! loop counter for nodes
integer(intKi) :: k ! loop counter for blades
Expand All @@ -3520,10 +3532,10 @@ subroutine Calculate_MeshOrientation_LiftingLine(p, u, m, thetaBladeNds, toeBlad

do j=1,p%NumBlNds
orientation = matmul( u%BladeMotion(k)%Orientation(:,:,j), transpose( m%orientationAnnulus(:,:,j,k) ) )
theta = EulerExtract( orientation )
m%Curve( j,k) = theta(2) ! TODO
if (present(thetaBladeNds)) thetaBladeNds(j,k) = -theta(3)
if (present(toeBladeNds )) toeBladeNds( j,k) = theta(1)
thetas = EulerExtract( orientation )
twist(j,k) = -thetas(3)
toe( j,k) = thetas(1)
cant( j,k) = thetas(2)
enddo
end do !k=blades

Expand Down Expand Up @@ -3563,15 +3575,16 @@ subroutine SetInputsForFVW(p, u, m, errStat, errMsg)
! 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, m%rotors(iR)%Toe, m%rotors(iR)%Cant, 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))
call TwistToeCant_FromLocalPolar(u(tIndx)%rotors(iR)%BladeMotion(k), m%rotors(iR)%orientationAnnulus(:,:,:,k), thetaBladeNds, m%rotors(iR)%Toe(:,k), m%rotors(iR)%Cant(:,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
call Calculate_MeshOrientation_LiftingLine (p%rotors(iR),u(tIndx)%rotors(iR), m%rotors(iR), thetaBladeNds, m%rotors(iR)%Toe, m%rotors(iR)%Cant, ErrStat=ErrStat2,ErrMsg=ErrMsg2) ! sets m%orientationAnnulus, m%Curve
else
call SetErrStat(ErrID_Fatal, 'Aero Projection Method not implemented' ,ErrStat, ErrMsg, RoutineName)
endif
Expand Down
4 changes: 2 additions & 2 deletions modules/aerodyn/src/AeroDyn_AllBldNdOuts_IO.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -672,7 +672,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, RotI
if (p_AD%Wake_Mod /= WakeMod_FVW) then
DO iB=1,nB
DO iNd=1,nNd
y%WriteOutput(iOut) = m%Curve(iNd,iB)*R2D
y%WriteOutput(iOut) = m%Cant(iNd,iB)*R2D
iOut = iOut + 1
END DO
END DO
Expand All @@ -691,7 +691,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, RotI
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
y%WriteOutput(iOut) = m%Toe(iNd,iB)*R2D
iOut = iOut + 1
END DO
END DO
Expand Down
7 changes: 3 additions & 4 deletions modules/aerodyn/src/AeroDyn_IO.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -140,6 +140,7 @@ SUBROUTINE Calc_WriteOutput( p, p_AD, u, RotInflow, x, m, m_AD, y, OtherState, x
endif


! Common outputs to all AeroDyn submodules
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%Wake_Mod /= WakeMod_FVW) then
Expand Down Expand Up @@ -236,7 +237,8 @@ subroutine Calc_WriteOutput_AD()
m%AllOuts( BNSTVy( beta,k) ) = tmp(2)
m%AllOuts( BNSTVz( beta,k) ) = tmp(3)

m%AllOuts( BNCurve(beta,k) ) = m%Curve(j,k)*R2D
m%AllOuts( BNCurve(beta,k) ) = m%Cant(j,k)*R2D
! TODO Toe

m%AllOuts( BNSigCr( beta,k) ) = m%SigmaCavitCrit(j,k)
m%AllOuts( BNSgCav( beta,k) ) = m%SigmaCavit(j,k)
Expand Down Expand Up @@ -422,8 +424,6 @@ subroutine Calc_WriteOutput_BEMT()
m%AllOuts( BNTheta(beta,k) ) = m%BEMT_u(indx)%theta(j,k)*R2D
m%AllOuts( BNPhi( beta,k) ) = m%BEMT_y%phi(j,k)*R2D

! m%AllOuts( BNCurve(beta,k) ) = m%Curve(j,k)*R2D

m%AllOuts( BNCpmin( beta,k) ) = m%BEMT_y%Cpmin(j,k)
! m%AllOuts( BNSigCr( beta,k) ) = m%SigmaCavitCrit(j,k)
! m%AllOuts( BNSgCav( beta,k) ) = m%SigmaCavit(j,k)
Expand Down Expand Up @@ -506,7 +506,6 @@ subroutine Calc_WriteOutput_FVW
m%AllOuts( BNAlpha(beta,k) ) = m_AD%FVW%W(iW)%BN_alpha(j)*R2D
m%AllOuts( BNTheta(beta,k) ) = m_AD%FVW%W(iW)%PitchAndTwist(j)*R2D
m%AllOuts( BNPhi( beta,k) ) = m_AD%FVW%W(iW)%BN_phi(j)*R2D
! m%AllOuts( BNCurve(beta,k) ) = m%Curve(j,k)*R2D ! TODO

m%AllOuts( BNCpmin(beta,k) ) = m_AD%FVW%W(iW)%BN_Cpmin(j)
m%AllOuts( BNCl( beta,k) ) = m_AD%FVW%W(iW)%BN_Cl(j)
Expand Down
3 changes: 2 additions & 1 deletion modules/aerodyn/src/AeroDyn_Registry.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -291,7 +291,8 @@ typedef ^ RotMiscVarType ReKi AllOuts {:} - - "An array holding the value of all
typedef ^ RotMiscVarType ReKi W_Twr {:} - - "relative wind speed normal to the tower at node j" m/s
typedef ^ RotMiscVarType ReKi X_Twr {:} - - "local x-component of force per unit length of the jth node in the tower" m/s
typedef ^ RotMiscVarType ReKi Y_Twr {:} - - "local y-component of force per unit length of the jth node in the tower" m/s
typedef ^ RotMiscVarType ReKi Curve {:}{:} - - "curvature angle, saved for possible output to file" rad
typedef ^ RotMiscVarType ReKi Cant {:}{:} - - "curvature angle, saved for possible output to file" rad
typedef ^ RotMiscVarType ReKi Toe {:}{:} - - "Toe angle, saved for possible output to file" rad
typedef ^ RotMiscVarType ReKi TwrClrnc {:}{:} - - "Distance between tower (including tower radius) and blade node (not including blade width), saved for possible output to file" m
typedef ^ RotMiscVarType ReKi X {:}{:} - - "normal force per unit length (normal to the plane, not chord) of the jth node in the kth blade" N/m
typedef ^ RotMiscVarType ReKi Y {:}{:} - - "tangential force per unit length (tangential to the plane, not chord) of the jth node in the kth blade" N/m
Expand Down
Loading

0 comments on commit 4430e82

Please sign in to comment.