Skip to content

Adding a simple nacelle drag model in AeroDyn #1888

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

Closed
mayankchetan wants to merge 4 commits into from
Closed

Adding a simple nacelle drag model in AeroDyn #1888

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
2b9332a
initial working nacelle drag
mayankchetan Nov 10, 2023
91d5692
new working approach
mayankchetan Nov 18, 2023
876f5cf
Cleaned up Nacelle Drag model, inc envision code in comments
mayankchetan Nov 22, 2023
139a8cb
few stylistic changes
mayankchetan Nov 22, 2023
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
Binary file modified docs/OtherSupporting/OutListParameters.xlsx
View file
Open in desktop
Binary file not shown.
153 changes: 150 additions & 3 deletions modules/aerodyn/src/AeroDyn.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1305,6 +1305,10 @@ subroutine SetParameters( InitInp, InputFileData, RotData, p, p_AD, ErrStat, Err
p%CavitCheck = InputFileData%CavitCheck
p%Buoyancy = InputFileData%Buoyancy

p%NacelleDrag = InputFileData%NacelleDrag
p%NacArea = RotData%NacArea
p%NacCd = RotData%NacCd
p%NacDragAC = RotData%NacDragAC

if (InitInp%Linearize .and. InputFileData%WakeMod == WakeMod_BEMT) then
p%FrozenWake = InputFileData%FrozenWake
Expand Down Expand Up @@ -1738,6 +1742,12 @@ subroutine AD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg,
call AD_CavtCrit(u, p, m, errStat2, errMsg2)
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)

! initialize nacelle mesh loads
do iR = 1,size(p%rotors)
y%rotors(iR)%NacelleLoad%Force = 0.0_ReKi
y%rotors(iR)%NacelleLoad%Moment = 0.0_ReKi
end do

! Calculate buoyant loads
do iR = 1,size(p%rotors)
if ( p%rotors(iR)%Buoyancy ) then
Expand All @@ -1747,6 +1757,15 @@ subroutine AD_CalcOutput( t, u, p, x, xd, z, OtherState, y, m, ErrStat, ErrMsg,
end if
end do

! Calculate nacelle drag loads
do iR = 1,size(p%rotors)
if ( p%rotors(iR)%NacelleDrag ) then
call computeNacelleDrag( u%rotors(iR), p%rotors(iR), m%rotors(iR), y%rotors(iR), ErrStat, ErrMsg )
call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
if (ErrStat >= AbortErrLev) return
end if
end do

!-------------------------------------------------------
! get values to output to file:
!-------------------------------------------------------
Expand Down Expand Up @@ -2352,9 +2371,9 @@ subroutine CalcBuoyantLoads( u, p, m, y, ErrStat, ErrMsg )
m%NacMB = NacMBtmp
end if

! Assign buoyant loads to nacelle mesh
y%NacelleLoad%Force(:,1) = NacFBtmp
y%NacelleLoad%Moment(:,1) = NacMBtmp
! Assign buoyant loads to nacelle mesh. Mesh might contain the nacelle drag force.
y%NacelleLoad%Force(:,1) = y%NacelleLoad%Force(:,1) + NacFBtmp
y%NacelleLoad%Moment(:,1) = y%NacelleLoad%Moment(:,1) + NacMBtmp

end subroutine CalcBuoyantLoads
!----------------------------------------------------------------------------------------------------------------------------------
Expand Down Expand Up @@ -3788,6 +3807,26 @@ SUBROUTINE ValidateInputData( InitInp, InputFileData, NumBl, ErrStat, ErrMsg )
call SetErrStat( ErrID_Fatal, 'DBEMT requires the continuous formulation with constant tau1 for linearization. Set DBEMT_Mod=3 or set WakeMod to 0 or 1.', ErrStat, ErrMsg, RoutineName )
end if
end if

if (InputFileData%NacelleDrag) then
call SetErrStat( ErrID_Fatal, 'Nacelle drag cannot currently be used for linearization. Set NacelleDrag = false.', ErrStat, ErrMsg, RoutineName )
end if
end if


!..................
! check for nacelle drag parameters
!..................

if (InputFileData%NacelleDrag) then
do iR = 1,size(NumBl)
if (any(InputFileData%rotors(iR)%NacArea < 0.0_ReKi)) then
call SetErrStat( ErrID_Fatal, 'Nacelle projected area should not be negative for drag model.', ErrStat, ErrMsg, RoutineName )
end if
if (any(InputFileData%rotors(iR)%NacCd < 0.0_ReKi)) then
call SetErrStat( ErrID_Fatal, 'Nacelle drag coefficient should not be negative for drag model.', ErrStat, ErrMsg, RoutineName )
end if
end do
end if

contains
Expand Down Expand Up @@ -7422,5 +7461,113 @@ subroutine AD_SetExternalWindPositions(u_AD, o_AD, PosXYZ, node, errStat, errMsg
enddo !j, wake points
end if
end subroutine AD_SetExternalWindPositions
!-------------------------------------------------------------------------------------------------------
!> This routine calculates nacelle drag loads on a turbine.
SUBROUTINE computeNacelleDrag( u, p, m, y, ErrStat, ErrMsg )

TYPE(RotInputType) , INTENT(IN ) :: u !< AD inputs - used for mesh node positions
TYPE(RotParameterType) , INTENT(IN ) :: p !< Parameters
TYPE(RotMiscVarType) , INTENT(INOUT) :: m !< Misc/optimization variables
TYPE(RotOutputType) , INTENT(INOUT) :: y !< Outputs computed at t
INTEGER(IntKi) , INTENT( OUT) :: ErrStat !< Error status of the operation
CHARACTER(*) , INTENT( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None
! Local Vars
REAL(ReKi) :: Nac_wnd(3) ! Wind velocity at nacelle mesh
REAL(ReKi) :: Nac_motion(3) ! Translational velocity of the nacelle
REAL(ReKi) :: Nac_Vrel(3) ! Relative velocity of the wind at nacelle
REAL(ReKi) :: Nac_Vrel_n(3) ! Relative wind velocity in the nacelle c.s
REAL(ReKi) :: Nac_Vrel_n_nom ! norm of realtive wind velocity
REAL(ReKi) :: f_drag_n(3) ! Drag force at nacelle mesh in nacelle c.s
REAL(ReKi) :: m_drag_n(3) ! Drag moment at nacelle mesh in nacelle c.s

! ! Vars for Envision function
! REAL(ReKi) :: totalAngle ! Angle between incoming wind direction and nacelle,
! REAL(ReKi) :: area ! Area of the nacelle projected in the wind direction
! REAL(ReKi) :: forceMag ! Drag force aligned with wind direction
! Real(ReKi) :: unitDiskVec(3) ! unit vector aligned at an angle of "totalAngle" from yawed rotor disk
! Real(ReKi) :: areaVec(3) ! Vec containing areas of yz, xz and xy faces of the nacelle
! REAL(ReKi) :: rho ! Air density
! REAL(ReKi) :: nacelleCD ! Nacelle Drag Coefficient
! REAL(ReKi) :: hubHeigthWindSpeed(3) ! hubHeigthWindSpeed(1), hubHeigthWindSpeed(2), and hubHeigthWindSpeed(3) and u, v, and w wind velocities at Hub height
! REAL(ReKi) :: nacelleDims(3) ! nacelleDims(1), nacelleDims(2), and nacelleDims(3) and Length(x), Width(y), and Height(z) of the Nacelle
! REAL(ReKi) :: yawAngle ! Current Yaw Bearing (Radians)
! REAL(ReKi) :: force(3) ! Forces in global c.s
! REAL(ReKi) :: moment(3) ! Moments in global c.s

ErrStat = ErrID_None
ErrMsg = ""

! Fetch nacelle inflow, and nacelle motion
Nac_wnd = u%InflowOnNacelle
Nac_motion = u%NacelleMotion%TranslationVel(:,1)

! Calculating the relative inflow velocity
Nac_Vrel = Nac_wnd - Nac_motion

! transforming from inertial to nacelle coordinate system (c.s)
Nac_Vrel_n = matmul(u%NacelleMotion%Orientation(:,:,1), Nac_Vrel)

! Calculating the norm of the relative inflow velocity
Nac_Vrel_n_nom = TwoNorm(Nac_Vrel_n)

! Find drag force in nacelle c.s, assuming projected area is also in nacelle c.s, acting at AC
f_drag_n(1:3) = 0.5 * p%AirDens * p%NacCd(1:3) * Nac_Vrel_n_nom * Nac_Vrel_n(1:3) * p%NacArea(1:3)

! moment affect due to offset between nacelle reference position and nacelle Drag AC
m_drag_n = CROSS_PRODUCT(p%NacDragAC, f_drag_n)

m%NacDragF = f_drag_n
m%NacDragM = m_drag_n

! Transfer to global c.s & adding the drag forces and moments to the Nacelle point mesh.
! Nacelle point mesh is currently defined at the tower top and will be at (0,0,TowerHt) for a landbased HAWT
y%NacelleLoad%Force(1:3,1) = y%NacelleLoad%Force(1:3,1) + matmul(transpose(u%NacelleMotion%Orientation(:,:,1)), f_drag_n)
y%NacelleLoad%Moment(1:3,1) = y%NacelleLoad%Moment(1:3,1) + matmul(transpose(u%NacelleMotion%Orientation(:,:,1)), m_drag_n)

! !! Code from Envision:
! ! Compatability with current subroutine.
! rho = p%AirDens
! nacelleCD = p%NacCd(1)
! hubHeigthWindSpeed = u%InflowOnNacelle
! nacelleDims = p%NacArea
! yawAngle = m%yaw

! totalAngle = atan2(hubHeigthWindSpeed(2),hubHeigthWindSpeed(1)) - yawAngle
! call MPi2Pi(totalAngle)
! unitDiskVec(1) = abs(cos(totalAngle))
! unitDiskVec(2) = abs(sin(totalAngle))
! unitDiskVec(3) = 0.0_ReKi ! Can be used if the tilt is to be included.

! areaVec(1) = nacelleDims(1)!nacelleDims(2)*nacelleDims(3) ! area as viewed from front (yz)
! areaVec(2) = nacelleDims(2)!nacelleDims(1)*nacelleDims(3) ! area as viewed from side (xz)
! areaVec(3) = nacelleDims(3)!nacelleDims(1)*nacelleDims(2) ! area as viewed from top (xy)

! ! total nacelle area projected into incoming wind direction
! area = dot_product(areaVec, unitDiskVec)

! ! Find drag force (in global X direction)
! forceMag = 0.5*rho*nacelleCD*(hubHeigthWindSpeed(1)**2 + hubHeigthWindSpeed(2)**2)*area

! ! Decompose along the nacelle length, width and height
! force = unitDiskVec*forceMag

! force(1) = sign(force(1),cos(totalAngle))
! force(2) = sign(force(2),sin(totalAngle))
! force(3) = 0.0_ReKi

! moment(1) = 0.0_ReKi
! moment(2) = 0.0_ReKi
! moment(3) = 0.0_ReKi

! Need to change from Inertial c.s to Naccelle c.s for outputs
! m%NacDragF = matmul(transpose(u%NacelleMotion%Orientation(:,:,1)),force)
! m%NacDragM = matmul(transpose(u%NacelleMotion%Orientation(:,:,1)),moment)

! y%NacelleLoad%Force(1:3,1) = y%NacelleLoad%Force(1:3,1) + force
! y%NacelleLoad%Moment(1:3,1) = y%NacelleLoad%Moment(1:3,1) + moment


END SUBROUTINE computeNacelleDrag

!----------------------------------------------------------------------------------------------------------------------------------
END MODULE AeroDyn
55 changes: 53 additions & 2 deletions modules/aerodyn/src/AeroDyn_IO.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -202,7 +202,7 @@ subroutine Calc_WriteOutput_AD()
m%AllOuts( HbMbz ) = tmpHubMB(3)
end if

! nacelle outputs
! nacelle buoyancy outputs
if ( p%Buoyancy ) then
tmp = matmul( u%NacelleMotion%Orientation(:,:,1) , m%NacFB )
m%AllOuts( NcFbx ) = tmp(1)
Expand All @@ -215,6 +215,24 @@ subroutine Calc_WriteOutput_AD()
m%AllOuts( NcMbz ) = tmp(3)
end if

! nacelle drag outputs
if ( p%NacelleDrag ) then
tmp = matmul( u%NacelleMotion%Orientation(:,:,1) , u%InflowOnNacelle )
m%AllOuts( NcVUndx ) = tmp(1)
m%AllOuts( NcVUndy ) = tmp(2)
m%AllOuts( NcVUndz ) = tmp(3)

tmp = m%NacDragF
m%AllOuts( NcFdx ) = tmp(1)
m%AllOuts( NcFdy ) = tmp(2)
m%AllOuts( NcFdz ) = tmp(3)

tmp = m%NacDragM
m%AllOuts( NcMdx ) = tmp(1)
m%AllOuts( NcMdy ) = tmp(2)
m%AllOuts( NcMdz ) = tmp(3)
end if

! blade outputs
do k=1,min(p%numBlades,AD_MaxBl_Out) ! limit this
do beta=1,p%NBlOuts
Expand Down Expand Up @@ -725,6 +743,9 @@ SUBROUTINE ParsePrimaryFileInfo( PriPath, InitInp, InputFile, RootName, NumBlade
if (Failed()) return
! Buoyancy - Include buoyancy effects? (flag)
call ParseVar( FileInfo_In, CurLine, "Buoyancy", InputFileData%Buoyancy, ErrStat2, ErrMsg2, UnEc )
if (Failed()) return
! NacelleDrag - Include Nacelle Drag effects? (flag)
call ParseVar( FileInfo_In, CurLine, "NacelleDrag", InputFileData%NacelleDrag, ErrStat2, ErrMsg2, UnEc )
if (Failed()) return
! CompAA - Flag to compute AeroAcoustics calculation [only used when WakeMod=1 or 2]
call ParseVar( FileInfo_In, CurLine, "CompAA", InputFileData%CompAA, ErrStat2, ErrMsg2, UnEc )
Expand Down Expand Up @@ -897,6 +918,16 @@ SUBROUTINE ParsePrimaryFileInfo( PriPath, InitInp, InputFile, RootName, NumBlade
! NacCenB - Nacelle center of buoyancy x,y,z direction offsets (m)
call ParseAry( FileInfo_In, CurLine, 'NacCenB', InputFileData%rotors(iR)%NacCenB, 3 , ErrStat2, ErrMsg2, UnEc )
if (Failed()) return

! NacArea - Projected area of the nacelle in X, Y, Z in the nacelle coordinate system (m^2)
call ParseAry( FileInfo_In, CurLine, "NacArea", InputFileData%rotors(iR)%NacArea, 3, ErrStat2, ErrMsg2, UnEc )
if (Failed()) return
! NacCd - Drag cefficient for the nacelle areas defied above (-)
call ParseAry( FileInfo_In, CurLine, "NacCd", InputFileData%rotors(iR)%NacCd, 3, ErrStat2, ErrMsg2, UnEc )
if (Failed()) return
! NacDragAC - Position of aerodynamic center of nacelle drag in nacelle coordinates (m)
call ParseAry( FileInfo_In, CurLine, "NacDragAC", InputFileData%rotors(iR)%NacDragAC, 3, ErrStat2, ErrMsg2, UnEc )
if (Failed()) return
end do

!====== Tail fin aerodynamics ========================================================================
Expand Down Expand Up @@ -1413,6 +1444,14 @@ SUBROUTINE AD_PrintSum( InputFileData, p, p_AD, u, y, ErrStat, ErrMsg )
end if
WRITE (UnSu,Ec_LgFrmt) p%Buoyancy, 'Buoyancy', 'Include buoyancy effects? '//TRIM(Msg)

! Nacelle Drag
if (p%NacelleDrag) then
Msg = 'Yes'
else
Msg = 'No'
end if
WRITE (UnSu,Ec_LgFrmt) p%NacelleDrag, 'NacelleDrag', 'Include NacelleDrag effects? '//TRIM(Msg)


if (p_AD%WakeMod/=WakeMod_none) then
WRITE (UnSu,'(A)') '====== Blade-Element/Momentum Theory Options ======================================================'
Expand Down Expand Up @@ -1693,7 +1732,19 @@ SUBROUTINE SetOutParam(OutList, p, p_AD, ErrStat, ErrMsg )
InvalidOutput( BNMbs(:,i) ) = .true.
end do
end if


if (.not. p%NacelleDrag) then ! Invalid Nacelle Drag loads
InvalidOutput( NcVUndx ) = .true.
InvalidOutput( NcVUndy ) = .true.
InvalidOutput( NcVUndz ) = .true.
InvalidOutput( NcFdx ) = .true.
InvalidOutput( NcFdy ) = .true.
InvalidOutput( NcFdz ) = .true.
InvalidOutput( NcMdx ) = .true.
InvalidOutput( NcMdy ) = .true.
InvalidOutput( NcMdz ) = .true.
end if

DO i = p%NTwOuts+1,9 ! Invalid tower nodes

InvalidOutput( TwNVUnd(:,i) ) = .true.
Expand Down
Loading