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.

Sign up for GitHub

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

[OLAF] Bug Fix: Twist Not Computed in OLAF cases when Polar Projection is used, resulting in wrong Fn, Ft outputs #2349

Merged
merged 1 commit into from
Aug 8, 2024
Merged

[OLAF] Bug Fix: Twist Not Computed in OLAF cases when Polar Projection is used, resulting in wrong Fn, Ft outputs #2349

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
99 changes: 56 additions & 43 deletions modules/aerodyn/src/AeroDyn.f90
View file
Open in desktop
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
41 changes: 10 additions & 31 deletions modules/aerodyn/src/AeroDyn_AllBldNdOuts_IO.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -669,41 +669,20 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, RotI
endif

CASE ( BldNd_Curve )
if (p_AD%Wake_Mod /= WakeMod_FVW) then
DO iB=1,nB
DO iNd=1,nNd
y%WriteOutput(iOut) = m%Curve(iNd,iB)*R2D
iOut = iOut + 1
END DO
END DO
else
DO iB=1,nB
iW = W2B(iB)
DO iNd=1,nNd
!NOT available in FVW yet
y%WriteOutput(iOut) = 0.0_ReKi
iOut = iOut + 1
END DO
DO iB=1,nB
DO iNd=1,nNd
y%WriteOutput(iOut) = m%Cant(iNd,iB)*R2D
iOut = iOut + 1
END DO
endif
END DO

CASE ( BldNd_Toe )
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
iOut = iOut + 1
END DO
END DO
else
DO iB=1,nB
iW = W2B(iB)
DO iNd=1,nNd
y%WriteOutput(iOut) = 0.0_ReKi
iOut = iOut + 1
END DO
DO iB=1,nB
DO iNd=1,nNd
y%WriteOutput(iOut) = m%Toe(iNd,iB)*R2D
iOut = iOut + 1
END DO
endif
END DO


! Unsteady lift force, drag force, pitching moment coefficients
Expand Down
6 changes: 2 additions & 4 deletions modules/aerodyn/src/AeroDyn_IO.f90
View file
Open in desktop
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,7 @@ 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

m%AllOuts( BNSigCr( beta,k) ) = m%SigmaCavitCrit(j,k)
m%AllOuts( BNSgCav( beta,k) ) = m%SigmaCavit(j,k)
Expand Down Expand Up @@ -422,8 +423,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 +505,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
Open in desktop
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