Skip to content

Commit

Permalink
Solve for induced velocity only in streamtubes where blades are locat…
Browse files Browse the repository at this point in the history 
…ed at each time step
  • Loading branch information
@hkross
hkross committed Oct 24, 2024
1 parent c95dd86 commit ec81bdf
Show file tree
Hide file tree
Showing 8 changed files with 253 additions and 325 deletions.
56 changes: 34 additions & 22 deletions modules/aerodyn/src/AeroDyn.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -459,8 +459,12 @@ subroutine AD_Init( InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut

elseif (p%Wake_Mod == WakeMod_DMST) then
do iR = 1, nRotors
call Init_DMSTmodule( InputFileData, InputFileData%rotors(iR), u%rotors(iR), m%rotors(iR)%DMST_u, p%rotors(iR), p, m%rotors(iR)%DMST_y, ErrStat2, ErrMsg2 )
call Init_DMSTmodule( InputFileData, InputFileData%rotors(iR), u%rotors(iR), m%rotors(iR)%DMST_u(1), p%rotors(iR), p, m%rotors(iR)%DMST_y, ErrStat2, ErrMsg2 )
if (Failed()) return;

call DMST_CopyInput( m%rotors(iR)%DMST_u(1), m%rotors(iR)%DMST_u(2), MESH_NEWCOPY, ErrStat2, ErrMsg2 )
call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )

enddo

else ! if (p%Wake_Mod == WakeMod_FVW) then
Expand Down Expand Up @@ -1841,9 +1845,17 @@ subroutine AD_UpdateStates( t, n, u, utimes, p, x, xd, z, OtherState, m, errStat
if (allocated(OtherState%WakeLocationPoints)) then
OtherState%WakeLocationPoints = m%FVW%r_wind
endif

elseif (p%Wake_Mod == WakeMod_DMST) then
do iR = 1,size(p%rotors)
call DMST_UpdateStates(m%rotors(iR)%DMST_u(:), errStat2, errMsg2)
if (Failed()) return
end do

! UA TODO
!call UA_UpdateState_Wrapper(p%AFI, n, p%FVW, x%FVW, xd%FVW, OtherState%FVW, m%FVW, ErrStat2, ErrMsg2)
! if (Failed()) return

endif

call Cleanup()
Expand Down Expand Up @@ -2177,7 +2189,7 @@ subroutine RotCalcOutput( t, u, RotInflow, p, p_AD, x, xd, z, OtherState, y, m,
elseif (p_AD%Wake_Mod == WakeMod_DMST) then
! Call the DMST module CalcOutput. Notice that the DMST outputs are purposely attached to AeroDyn's MiscVar structure to
! avoid issues with the coupling code
call DMST_CalcOutput( m%DMST_u, p%DMST, p_AD%AFI, m%DMST_y, ErrStat2, ErrMsg2 )
call DMST_CalcOutput( m%DMST_u(indx), p%DMST, p_AD%AFI, m%DMST_y, ErrStat2, ErrMsg2 )
call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
call SetOutputsFromDMST( u, p, p_AD, m, y, ErrStat2, ErrMsg2 )
call SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
Expand Down Expand Up @@ -2879,7 +2891,7 @@ subroutine SetInputs(t, p, p_AD, u, RotInflow, m, indx, errStat, errMsg)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
elseif (p_AD%Wake_Mod == WakeMod_DMST) then
! This needs to extract the inputs from the AD data types (mesh) and massage them for the DMST module
call SetInputsForDMST(p_AD, p, u, RotInflow, m, errStat2, errMsg2)
call SetInputsForDMST(p, p_AD, u, RotInflow, m, indx, errStat2, errMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
endif

Expand Down Expand Up @@ -3652,13 +3664,14 @@ subroutine Calculate_MeshOrientation_LiftingLine(p, u, m, twist, toe, cant, ErrS
end subroutine Calculate_MeshOrientation_LiftingLine
!----------------------------------------------------------------------------------------------------------------------------------
!> This subroutine sets m%DMST_u.
subroutine SetInputsForDMST(p_AD, p, u, RotInflow, m, errStat, errMsg)
subroutine SetInputsForDMST(p, p_AD, u, RotInflow, m, indx, errStat, errMsg)

type(AD_ParameterType), intent(in ) :: p_AD !< AD parameters
type(RotParameterType), intent(in ) :: p !< AD parameters
type(AD_ParameterType), intent(in ) :: p_AD !< AD parameters
type(RotInputType), intent(in ) :: u !< AD inputs at time
type(RotInflowType), intent(in ) :: RotInflow !< Rotor inflow at time
type(RotMiscVarType), intent(inout) :: m !< Misc/optimization variables
integer, intent(in ) :: indx !< index into m%DMST_u array; must be 1 or 2 (but not checked here)
integer(IntKi), intent( out) :: ErrStat !< Error status of the operation
character(*), intent( out) :: ErrMsg !< Error message if ErrStat /= ErrID_None

Expand All @@ -3680,18 +3693,14 @@ subroutine SetInputsForDMST(p_AD, p, u, RotInflow, m, errStat, errMsg)

call Calculate_MeshOrientation_LiftingLine( p, u, m, thetaBladeNds, m%Toe, m%Cant, ErrStat=ErrStat, ErrMsg=ErrMsg )
if (ErrStat >= AbortErrLev) return
m%DMST_u%PitchAndTwist = thetaBladeNds
m%DMST_u(indx)%PitchAndTwist = thetaBladeNds
if (allocated(thetaBladeNds)) deallocate(thetaBladeNds)

! Free-stream velocity, m/s
do j = 1,p%NumBlNds
do i = 1,p%DMST%Nst
m%DMST_u%Vinf(:,i,j) = m%DisturbedInflow(:,j,1)
end do
end do
m%DMST_u(indx)%Vinf = m%DisturbedInflow

! Rotor angular velocity, rad/s
m%DMST_u%omega = dot_product( u%HubMotion%Orientation(1,1:3,1), u%HubMotion%RotationVel(:,1) )
m%DMST_u(indx)%omega = dot_product( u%HubMotion%Orientation(1,1:3,1), u%HubMotion%RotationVel(:,1) )

! Azimuthal location of each blade, rad
theta_b_tmp = EulerExtract( matmul(u%HubMotion%Orientation(:,:,1), transpose(u%HubMotion%RefOrientation(:,:,1))) )
Expand All @@ -3705,8 +3714,8 @@ subroutine SetInputsForDMST(p_AD, p, u, RotInflow, m, errStat, errMsg)
! Range of azimuth angles within each streamtube
do j = 1,p%NumBlNds
do i = 1,2*p%DMST%Nst
theta_st_r(1,i,j) = p%DMST%theta_st(i,j) - p%DMST%dTheta(j)/2.0_ReKi
theta_st_r(2,i,j) = p%DMST%theta_st(i,j) + p%DMST%dTheta(j)/2.0_ReKi
theta_st_r(1,i,j) = p%DMST%theta_st(i) - p%DMST%dTheta(j)/2.0_ReKi
theta_st_r(2,i,j) = p%DMST%theta_st(i) + p%DMST%dTheta(j)/2.0_ReKi
end do
end do

Expand All @@ -3715,13 +3724,13 @@ subroutine SetInputsForDMST(p_AD, p, u, RotInflow, m, errStat, errMsg)
do j = 1,p%NumBlNds
do i = 1,2*p%DMST%Nst
if ( theta_b(k) >= theta_st_r(1,i,j) .and. theta_b(k) < theta_st_r(2,i,j) ) then
m%DMST_u%blade_st(j,k) = i
m%DMST_u(indx)%blade_st(j,k) = i
end if
end do
end do
end do

m%DMST_u%UserProp = u%UserProp
m%DMST_u(indx)%UserProp = u%UserProp

end subroutine SetInputsForDMST
!----------------------------------------------------------------------------------------------------------------------------------
Expand Down Expand Up @@ -4103,6 +4112,7 @@ subroutine SetOutputsFromDMST(u, p, p_AD, m, y, ErrStat, ErrMsg)
real(reki) :: moment(3)
real(reki) :: q
real(ReKi) :: cp, sp ! Cosine, sine of phi
integer, parameter :: indx = 1

! Local variables for readability
real(ReKi) :: Vind(3)
Expand Down Expand Up @@ -4133,7 +4143,7 @@ subroutine SetOutputsFromDMST(u, p, p_AD, m, y, ErrStat, ErrMsg)
Vind = m%DMST_y%Vind(1:3,j,k)
Vstr = u%BladeMotion(k)%TranslationVel(1:3,j)
Vwnd = m%DisturbedInflow(1:3,j,k) ! contains tower shadow
theta = m%DMST_u%PitchAndTwist(j,k)
theta = m%DMST_u(indx)%PitchAndTwist(j,k)
call LL_AeroOuts( m%orientationAnnulus(1:3,1:3,j,k), u%BladeMotion(k)%Orientation(1:3,1:3,j), & ! inputs
theta, Vstr(1:3), Vind(1:3), Vwnd(1:3), p%KinVisc, p%DMST%chord(j,k), & ! inputs
AxInd, TanInd, Vrel, phi, alpha, Re, UrelWind_s(1:3), ErrStat2, ErrMsg2 ) ! outputs
Expand Down Expand Up @@ -5060,7 +5070,7 @@ SUBROUTINE Init_DMSTmodule( InputFileData, RotInputFileData, u_AD, u, p, p_AD, y
integer(intKi) :: k ! Blade index
real(ReKi) :: y_hat_disk(3)
real(ReKi) :: z_hat_disk(3)
real(ReKi) :: dr_gl(3)
real(ReKi) :: dr_gl(3,p%NumBlNds,p%NumBlades)
integer(IntKi) :: ErrStat2
character(ErrMsgLen) :: ErrMsg2
character(*), parameter :: RoutineName = 'Init_DMSTmodule'
Expand All @@ -5082,7 +5092,7 @@ SUBROUTINE Init_DMSTmodule( InputFileData, RotInputFileData, u_AD, u, p, p_AD, y

call AllocAry(InitInp%chord, InitInp%numBladeNodes,InitInp%numBlades,'chord', ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
call AllocAry(InitInp%AFindx,InitInp%numBladeNodes,InitInp%numBlades,'AFindx',ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
call AllocAry(InitInp%radius,InitInp%numBladeNodes,'radius',ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)
call AllocAry(InitInp%radius,InitInp%numBladeNodes,InitInp%numBlades,'radius',ErrStat2,ErrMsg2); call SetErrStat(ErrStat2,ErrMsg2,ErrStat,ErrMsg,RoutineName)

if ( ErrStat >= AbortErrLev ) then
call Cleanup()
Expand All @@ -5098,9 +5108,11 @@ SUBROUTINE Init_DMSTmodule( InputFileData, RotInputFileData, u_AD, u, p, p_AD, y

y_hat_disk = u_AD%HubMotion%Orientation(2,:,1)
z_hat_disk = u_AD%HubMotion%Orientation(3,:,1)
do j=1,p%NumBlNds
dr_gl = u_AD%BladeMotion(1)%Position(:,j) - u_AD%HubMotion%Position(:,1) ! vector hub center to node j in global coord
InitInp%radius(j) = sqrt( dot_product(dr_gl, y_hat_disk)**2 + dot_product(dr_gl, z_hat_disk)**2 )
do k=1,p%NumBlades
do j=1,p%NumBlNds
dr_gl(:,j,k) = u_AD%BladeMotion(k)%Position(:,j) - u_AD%HubMotion%Position(:,1) ! vector hub center to node j in global coord
InitInp%radius(j,k) = sqrt( dot_product(dr_gl(:,j,k), y_hat_disk)**2 + dot_product(dr_gl(:,j,k), z_hat_disk)**2 )
end do
end do

if (ErrStat >= AbortErrLev) then
Expand Down
2 changes: 1 addition & 1 deletion modules/aerodyn/src/AeroDyn_AllBldNdOuts_IO.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -644,7 +644,7 @@ SUBROUTINE Calc_WriteAllBldNdOutput( p, p_AD, u, m, m_AD, x, y, OtherState, RotI
else
DO iB=1,nB
DO iNdL=1,nNd; iNd=Nd(iNdL);
y%WriteOutput(iOut) = m_AD%rotors(iRot)%DMST_u%PitchAndTwist(iNd,iB)*R2D
y%WriteOutput(iOut) = m_AD%rotors(iRot)%DMST_u(indx)%PitchAndTwist(iNd,iB)*R2D
iOut = iOut + 1
END DO
END DO
Expand Down
2 changes: 1 addition & 1 deletion modules/aerodyn/src/AeroDyn_IO.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -561,7 +561,7 @@ subroutine Calc_WriteOutput_LL
if ( p_AD%Wake_Mod == WakeMod_FVW ) then
m%AllOuts( BNTheta(beta,k) ) = m_AD%FVW%W(iW)%PitchAndTwist(j)*R2D
else
m%AllOuts( BNTheta(beta,k) ) = m_AD%rotors(iRot)%DMST_u%PitchAndTwist(j,k)*R2D
m%AllOuts( BNTheta(beta,k) ) = m_AD%rotors(iRot)%DMST_u(indx)%PitchAndTwist(j,k)*R2D
endif
m%AllOuts( BNPhi( beta,k) ) = m_AD%rotors(iRot)%blds(k)%BN_phi(j)*R2D

Expand Down
2 changes: 1 addition & 1 deletion modules/aerodyn/src/AeroDyn_Registry.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -310,7 +310,7 @@ typedef ^ RotMiscVarType BEMT_MiscVarType BEMT - - - "MiscVars from the BEMT mod
typedef ^ RotMiscVarType BEMT_OutputType BEMT_y - - - "Outputs from the BEMT module" -
typedef ^ RotMiscVarType BEMT_InputType BEMT_u 2 - - "Inputs to the BEMT module" -
typedef ^ RotMiscVarType DMST_OutputType DMST_y - - - "Outputs from the DMST module" -
typedef ^ RotMiscVarType DMST_InputType DMST_u - - - "Inputs to the DMST module" -
typedef ^ RotMiscVarType DMST_InputType DMST_u 2 - - "Inputs to the DMST module" -
typedef ^ RotMiscVarType AA_MiscVarType AA - - - "MiscVars from the AA module" -
typedef ^ RotMiscVarType AA_OutputType AA_y - - - "Outputs from the AA module" -
typedef ^ RotMiscVarType AA_InputType AA_u - - - "Inputs to the AA module" -
Expand Down
32 changes: 24 additions & 8 deletions modules/aerodyn/src/AeroDyn_Types.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -351,7 +351,7 @@ MODULE AeroDyn_Types
TYPE(BEMT_OutputType) :: BEMT_y !< Outputs from the BEMT module [-]
TYPE(BEMT_InputType) , DIMENSION(1:2) :: BEMT_u !< Inputs to the BEMT module [-]
TYPE(DMST_OutputType) :: DMST_y !< Outputs from the DMST module [-]
TYPE(DMST_InputType) :: DMST_u !< Inputs to the DMST module [-]
TYPE(DMST_InputType) , DIMENSION(1:2) :: DMST_u !< Inputs to the DMST module [-]
TYPE(AA_MiscVarType) :: AA !< MiscVars from the AA module [-]
TYPE(AA_OutputType) :: AA_y !< Outputs from the AA module [-]
TYPE(AA_InputType) :: AA_u !< Inputs to the AA module [-]
Expand Down Expand Up @@ -3468,9 +3468,13 @@ subroutine AD_CopyRotMiscVarType(SrcRotMiscVarTypeData, DstRotMiscVarTypeData, C
call DMST_CopyOutput(SrcRotMiscVarTypeData%DMST_y, DstRotMiscVarTypeData%DMST_y, CtrlCode, ErrStat2, ErrMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
if (ErrStat >= AbortErrLev) return
call DMST_CopyInput(SrcRotMiscVarTypeData%DMST_u, DstRotMiscVarTypeData%DMST_u, CtrlCode, ErrStat2, ErrMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
if (ErrStat >= AbortErrLev) return
LB(1:1) = lbound(SrcRotMiscVarTypeData%DMST_u, kind=B8Ki)
UB(1:1) = ubound(SrcRotMiscVarTypeData%DMST_u, kind=B8Ki)
do i1 = LB(1), UB(1)
call DMST_CopyInput(SrcRotMiscVarTypeData%DMST_u(i1), DstRotMiscVarTypeData%DMST_u(i1), CtrlCode, ErrStat2, ErrMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
if (ErrStat >= AbortErrLev) return
end do
call AA_CopyMisc(SrcRotMiscVarTypeData%AA, DstRotMiscVarTypeData%AA, CtrlCode, ErrStat2, ErrMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
if (ErrStat >= AbortErrLev) return
Expand Down Expand Up @@ -4034,8 +4038,12 @@ subroutine AD_DestroyRotMiscVarType(RotMiscVarTypeData, ErrStat, ErrMsg)
end do
call DMST_DestroyOutput(RotMiscVarTypeData%DMST_y, ErrStat2, ErrMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
call DMST_DestroyInput(RotMiscVarTypeData%DMST_u, ErrStat2, ErrMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
LB(1:1) = lbound(RotMiscVarTypeData%DMST_u, kind=B8Ki)
UB(1:1) = ubound(RotMiscVarTypeData%DMST_u, kind=B8Ki)
do i1 = LB(1), UB(1)
call DMST_DestroyInput(RotMiscVarTypeData%DMST_u(i1), ErrStat2, ErrMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
end do
call AA_DestroyMisc(RotMiscVarTypeData%AA, ErrStat2, ErrMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
call AA_DestroyOutput(RotMiscVarTypeData%AA_y, ErrStat2, ErrMsg2)
Expand Down Expand Up @@ -4229,7 +4237,11 @@ subroutine AD_PackRotMiscVarType(RF, Indata)
call BEMT_PackInput(RF, InData%BEMT_u(i1))
end do
call DMST_PackOutput(RF, InData%DMST_y)
call DMST_PackInput(RF, InData%DMST_u)
LB(1:1) = lbound(InData%DMST_u, kind=B8Ki)
UB(1:1) = ubound(InData%DMST_u, kind=B8Ki)
do i1 = LB(1), UB(1)
call DMST_PackInput(RF, InData%DMST_u(i1))
end do
call AA_PackMisc(RF, InData%AA)
call AA_PackOutput(RF, InData%AA_y)
call AA_PackInput(RF, InData%AA_u)
Expand Down Expand Up @@ -4373,7 +4385,11 @@ subroutine AD_UnPackRotMiscVarType(RF, OutData)
call BEMT_UnpackInput(RF, OutData%BEMT_u(i1)) ! BEMT_u
end do
call DMST_UnpackOutput(RF, OutData%DMST_y) ! DMST_y
call DMST_UnpackInput(RF, OutData%DMST_u) ! DMST_u
LB(1:1) = lbound(OutData%DMST_u, kind=B8Ki)
UB(1:1) = ubound(OutData%DMST_u, kind=B8Ki)
do i1 = LB(1), UB(1)
call DMST_UnpackInput(RF, OutData%DMST_u(i1)) ! DMST_u
end do
call AA_UnpackMisc(RF, OutData%AA) ! AA
call AA_UnpackOutput(RF, OutData%AA_y) ! AA_y
call AA_UnpackInput(RF, OutData%AA_u) ! AA_u
Expand Down
Loading

0 comments on commit ec81bdf

Please sign in to comment.