Skip to content

beamdyn and orcaflex updates #38

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
ghost wants to merge 5 commits into from
Closed

beamdyn and orcaflex updates #38

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
fd8e84e
bjonkman/f/Envision-BeamDyn: the BD input file specified in the drive…
rafmudaf Jul 27, 2017
062941b
bjonkman/f/Envision-BeamDyn: rotational velocity vector is now double…
rafmudaf Jul 27, 2017
d76d4b4
use R8Ki in BeamDyn_Driver for support with single precision
rafmudaf Jul 27, 2017
6fc9086
bjonkman/f/Envision-BeamDyn: BD: updated BD_CrvExtractCrv for more ro…
rafmudaf Jul 27, 2017
f65ba70
use ReKi in OrcaDriver_Subs for support with single precision
rafmudaf Jul 28, 2017
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
100 changes: 46 additions & 54 deletions modules-local/beamdyn/src/BeamDyn_Subs.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -339,15 +339,10 @@ SUBROUTINE BD_CrvExtractCrv(R, cc, ErrStat, ErrMsg)
CHARACTER(*), INTENT( OUT) :: ErrMsg !< Error message if ErrStat /= ErrID_None

!Local variables
REAL(BDKi) :: pivot(4) ! Trace of the rotation matrix and diagonal elements
REAL(BDKi) :: sm0
REAL(BDKi) :: sm1
REAL(BDKi) :: sm2
REAL(BDKi) :: sm3
REAL(BDKi) :: em
REAL(BDKi) :: temp
integer :: ipivot(1)
REAL(BDKi) :: sm(0:3)
REAL(BDKi) :: em
REAL(BDKi) :: Rr(3,3) ! (possibly) corrected rotation matrix
INTEGER :: i ! case indicator

INTEGER(IntKi) :: ErrStat2 ! Temporary Error status
CHARACTER(ErrMsgLen) :: ErrMsg2 ! Temporary Error message
Expand Down Expand Up @@ -376,54 +371,51 @@ SUBROUTINE BD_CrvExtractCrv(R, cc, ErrStat, ErrMsg)
! mjs--determine whether R is a valid rotation matrix and correct it if not
call BD_CheckRotMat(R, Rr, ErrStat2, ErrMsg2)
CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
if (ErrStat >= AbortErrLev) return
if (ErrStat >= AbortErrLev) return

! mjs--find max value of T := Tr(Rr) and diagonal elements of Rr
pivot = (/ Rr(1, 1), Rr(2, 2), Rr(3, 3), Rr(1,1) + Rr(2,2) + Rr(3,3) /)
i = maxloc(pivot, 1)

!> - Condition 1: \f$ \underline{\underline{R(3,3)}} > T \f$
!! This implies that \f$ c_3^2 > c_0^2 \f$
IF (i == 3) THEN
sm0 = Rr(2,1) - Rr(1,2) ! 4 c_0 c_3 t_{r0}
sm1 = Rr(1,3) + Rr(3,1) ! 4 c_1 c_3 t_{r0}
sm2 = Rr(2,3) + Rr(3,2) ! 4 c_2 c_3 t_{r0}
sm3 = 1.0_BDKi - Rr(1,1) - Rr(2,2) + Rr(3,3) ! 4 c_3^2 t_{r0}
temp = SIGN( 2.0_BDKi*SQRT(ABS(sm3)), sm0 )

!> - Condition 2: \f$ \underline{\underline{R(2,2)}} > T \f$
!! This implies that \f$ c_2^2 > c_0^2 \f$
ELSEIF (i == 2) THEN
sm0 = Rr(1,3) - Rr(3,1) ! 4 c_0 c_2 t_{r0}
sm1 = Rr(1,2) + Rr(2,1) ! 4 c_1 c_2 t_{r0}
sm2 = 1.0_BDKi - Rr(1,1) + Rr(2,2) - Rr(3,3) ! 4 c_2^2 t_{r0}
sm3 = Rr(2,3) + Rr(3,2) ! 4 c_2 c_3 t_{r0}
temp = SIGN( 2.0_BDKi*SQRT(ABS(sm2)), sm0 )

!> - Condition 3: \f$ \underline{\underline{R(1,1)}} > T \f$
!! This implies that \f$ c_1^2 > c_0^2 \f$
ELSEIF (i == 1) THEN
sm0 = Rr(3,2) - Rr(2,3) ! 4 c_0 c_1 t_{r0}
sm1 = 1.0_BDKi + Rr(1,1) - Rr(2,2) - Rr(3,3) ! 4 c_1^2 t_{r0}
sm2 = Rr(1,2) + Rr(2,1) ! 4 c_1 c_2 t_{r0}
sm3 = Rr(1,3) + Rr(3,1) ! 4 c_1 c_3 t_{r0}
temp = SIGN( 2.0_BDKi*SQRT(ABS(sm1)), sm0 )

!> - Condition 4: all diagonal terms are less than the trace
ELSE
! Note that 0 <= c_0 <= 2 since 0 <= dot_product(cc,cc) <= 16

! Let's first find out which parameter is largest so we can divide by it later
sm(0) = 1.0_BDKi + Rr(1,1) + Rr(2,2) + Rr(3,3) ! 4 c_0 c_0 t_{r0}
sm(1) = 1.0_BDKi + Rr(1,1) - Rr(2,2) - Rr(3,3) ! 4 c_1 c_1 t_{r0}
sm(2) = 1.0_BDKi - Rr(1,1) + Rr(2,2) - Rr(3,3) ! 4 c_2 c_2 t_{r0}
sm(3) = 1.0_BDKi - Rr(1,1) - Rr(2,2) + Rr(3,3) ! 4 c_3 c_3 t_{r0}
ipivot = maxloc(sm) - 1 ! our sm array starts at 0, so we need to subtract 1 here to get the correct index

select case (ipivot(1))
case (3)
!! We need the condition that c_3 is not zero.
sm(0) = Rr(2,1) - Rr(1,2) ! 4 c_0 c_3 t_{r0}
sm(1) = Rr(1,3) + Rr(3,1) ! 4 c_1 c_3 t_{r0}
sm(2) = Rr(2,3) + Rr(3,2) ! 4 c_2 c_3 t_{r0}
!sm(3) = 1.0_BDKi - Rr(1,1) - Rr(2,2) + Rr(3,3) ! 4 c_3 c_3 t_{r0}

case (2)
!! We need the condition that c_2 is not zero.
sm(0) = Rr(1,3) - Rr(3,1) ! 4 c_0 c_2 t_{r0}
sm(1) = Rr(1,2) + Rr(2,1) ! 4 c_1 c_2 t_{r0}
!sm(2) = 1.0_BDKi - Rr(1,1) + Rr(2,2) - Rr(3,3) ! 4 c_2 c_2 t_{r0}
sm(3) = Rr(2,3) + Rr(3,2) ! 4 c_3 c_2 t_{r0}

case (1)
!! We need the condition that c_1 is not zero.
sm(0) = Rr(3,2) - Rr(2,3) ! 4 c_0 c_1 t_{r0}
!sm(1) = 1.0_BDKi + Rr(1,1) - Rr(2,2) - Rr(3,3) ! 4 c_1 c_1 t_{r0}
sm(2) = Rr(1,2) + Rr(2,1) ! 4 c_2 c_1 t_{r0}
sm(3) = Rr(1,3) + Rr(3,1) ! 4 c_3 c_1 t_{r0}

case (0)
!! We need the condition that c_0 is not zero.
!ipivot = 0
sm0 = 1.0_BDKi + Rr(1,1) + Rr(2,2) + Rr(3,3) ! 4 c_0^2 t_{r0}
sm1 = Rr(3,2) - Rr(2,3) ! -4 c_0 c_1 t_{r0}
sm2 = Rr(1,3) - Rr(3,1) ! -4 c_0 c_2 t_{r0}
sm3 = Rr(2,1) - Rr(1,2) ! -4 c_0 c_3 t_{r0}
temp = SIGN( 2.0_BDKi*SQRT(ABS(sm0)), sm0 )
ENDIF

em = sm0 + temp
em = 4.0_BDKi/em
cc(1) = em*sm1
cc(2) = em*sm2
cc(3) = em*sm3
!sm(0) = 1.0_BDKi + Rr(1,1) + Rr(2,2) + Rr(3,3) ! 4 c_0 c_0 t_{r0}
sm(1) = Rr(3,2) - Rr(2,3) ! 4 c_1 c_0 t_{r0}
sm(2) = Rr(1,3) - Rr(3,1) ! 4 c_2 c_0 t_{r0}
sm(3) = Rr(2,1) - Rr(1,2) ! 4 c_3 c_0 t_{r0}
end select

em = sm(0) + SIGN( 2.0_BDKi*SQRT(sm(ipivot(1))), sm(0) )
em = 4.0_BDKi/em ! 1 / ( 4 t_{r0} c_{ipivot} ), assuming 0 <= c_0 < 4 and c_{ipivot} > 0
cc = em*sm(1:3)

END SUBROUTINE BD_CrvExtractCrv
!------------------------------------------------------------------------------
Expand Down
53 changes: 27 additions & 26 deletions modules-local/beamdyn/src/Driver_Beam_Subs.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,26 +27,22 @@ module BeamDyn_driver_subs
TYPE , PUBLIC :: BD_DriverInternalType
REAL(ReKi) , DIMENSION(1:6) :: DistrLoad !< Constant distributed load along beam axis, 3 forces and 3 moments [-]
REAL(ReKi) , DIMENSION(1:6) :: TipLoad !< Constant point load applied at tip, 3 forces and 3 moments [-]
INTEGER(IntKi) :: NumPointLoads !< Number of constant point loads applied along beam axis, 3 forces and 3 moments [-]
REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE :: MultiPointLoad !< Constant point loads applied along beam axis (index 1= Relative position along blade span; indices 2-7 = Fx, Fy, Fz, Mx, My, Mz) [-]
! INTEGER(IntKi) :: NumPointLoads !< Number of multi-point loads in the driver input file
! REAL(DbKi) :: MultiPointLoads{:}{:} !< The array of multipoint loads Index 1: [1, NumPointLoads];
!< Index 2: [1,7] (index of Loads 1 = Relative position along blade span;
!! 2-7 = Fx, Fy, Fz, Mx, My, Mz )

TYPE(MeshType) :: mplMotion ! Mesh for blade motion at multipoint loads locations
TYPE(MeshType) :: mplLoads ! Mesh for multipoint loads
TYPE(MeshMapType) :: Map_BldMotion_to_mplMotion
TYPE(MeshMapType) :: Map_mplLoads_to_PointLoad
TYPE(MeshType) :: y_BldMotion_at_u_point ! Intermediate mesh to transfer motion from output mesh to input mesh
TYPE(MeshMapType) :: Map_y_BldMotion_to_u_point

TYPE(MeshType) :: RotationCenter
TYPE(MeshMapType) :: Map_RotationCenter_to_RootMotion
INTEGER(IntKi) :: NumPointLoads !< Number of constant point loads applied along beam axis, 3 forces and 3 moments, from the driver input file [-]
REAL(BDKi) , DIMENSION(:,:), ALLOCATABLE :: MultiPointLoad !< Constant point loads applied along beam axis, size (NumPointLoads,7); (dimension 2: index 1=Relative position along blade span; indices 2-7 = Fx, Fy, Fz, Mx, My, Mz) [-]

TYPE(MeshType) :: mplMotion ! Mesh for blade motion at multipoint loads locations
TYPE(MeshType) :: mplLoads ! Mesh for multipoint loads
TYPE(MeshMapType) :: Map_BldMotion_to_mplMotion
TYPE(MeshMapType) :: Map_mplLoads_to_PointLoad
TYPE(MeshType) :: y_BldMotion_at_u_point ! Intermediate mesh to transfer motion from output mesh to input mesh
TYPE(MeshMapType) :: Map_y_BldMotion_to_u_point

REAL(DbKi) :: t_initial
REAL(DbKi) :: t_final
REAL(ReKi) :: w ! magnitude of rotational velocity vector
TYPE(MeshType) :: RotationCenter
TYPE(MeshMapType) :: Map_RotationCenter_to_RootMotion

REAL(DbKi) :: t_initial
REAL(DbKi) :: t_final
REAL(R8Ki) :: w ! magnitude of rotational velocity vector

END TYPE

Expand Down Expand Up @@ -83,6 +79,7 @@ SUBROUTINE BD_ReadDvrFile(DvrInputFile,dt,InitInputData,DvrData,&
INTEGER(IntKi) :: UnEc

CHARACTER(1024) :: FTitle ! "File Title": the 2nd line of the input file, which contains a description of its contents
CHARACTER(1024) :: PriPath ! Path name of the primary file

INTEGER(IntKi) :: i
!------------------------------------------------------------------------------------
Expand All @@ -98,6 +95,8 @@ SUBROUTINE BD_ReadDvrFile(DvrInputFile,dt,InitInputData,DvrData,&
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
IF (ErrStat >= AbortErrLev) RETURN

CALL GetPath( DvrInputFile, PriPath ) ! Input files will be relative to the path where the primary input file is located.

!-------------------------- HEADER ---------------------------------------------
CALL ReadCom(UnIn,DvrInputFile,'File Header: Module Version (line 1)',ErrStat2,ErrMsg2,UnEc)
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
Expand Down Expand Up @@ -270,6 +269,8 @@ SUBROUTINE BD_ReadDvrFile(DvrInputFile,dt,InitInputData,DvrData,&
!---------------------- BEAM SECTIONAL PARAMETER ----------------------------------------
CALL ReadVar ( UnIn, DvrInputFile, InitInputData%InputFile, 'InputFile', 'Name of the primary input file', ErrStat2,ErrMsg2, UnEc )
CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
IF ( PathIsRelative( InitInputData%InputFile ) ) InitInputData%InputFile = TRIM(PriPath)//TRIM(InitInputData%InputFile)

if (ErrStat >= AbortErrLev) then
call cleanup()
return
Expand Down Expand Up @@ -586,20 +587,20 @@ SUBROUTINE Init_RotationCenterMesh(DvrData, InitInputData, RootMotionMesh, ErrSt

DvrData%w = TwoNorm( InitInputData%RootVel(4:6) )

if (EqualRealNos(DvrData%w,0.0_ReKi)) then
DvrData%w = 0.0_ReKi
if (EqualRealNos(DvrData%w,0.0_R8Ki)) then
DvrData%w = 0.0_R8Ki
! the beam is not rotating, so pick an orientation
call eye(orientation, ErrStat2, ErrMsg2)
else
z_hat = InitInputData%RootVel(4:6) / DvrData%w

if ( EqualRealNos( z_hat(3), 1.0_ReKi ) ) then
if ( EqualRealNos( z_hat(3), 1.0_R8Ki ) ) then
call eye(orientation, ErrStat2, ErrMsg2)
elseif ( EqualRealNos( z_hat(3), -1.0_ReKi ) ) then
elseif ( EqualRealNos( z_hat(3), -1.0_R8Ki ) ) then
orientation = 0.0_ReKi
orientation(1,1) = -1.0_ReKi
orientation(2,2) = 1.0_ReKi
orientation(3,3) = -1.0_ReKi
orientation(1,1) = -1.0_R8Ki
orientation(2,2) = 1.0_R8Ki
orientation(3,3) = -1.0_R8Ki
else

Z_unit = (/0, 0, 1/)
Expand Down
43 changes: 21 additions & 22 deletions modules-local/orcaflex-interface/src/OrcaDriver_Subs.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1907,12 +1907,13 @@ SUBROUTINE PointsForce_OutputWrite(ProgInfo, OutUnit, OutFileName, InputFileName

! Temporary local variables
INTEGER(IntKi) :: ErrStatTmp !< Temporary variable for the status of error message
CHARACTER(*), PARAMETER :: RoutineName = 'PointsForce_OutputWrite'
CHARACTER(*), PARAMETER :: RoutineName = 'PointsForce_OutputWrite'
CHARACTER(2048) :: ErrMsgTmp !< Temporary variable for the error message
INTEGER(IntKi) :: LenErrMsgTmp !< Length of ErrMsgTmp (for getting WindGrid info)
REAL(ReKi) :: rotdisp(3) !< Rotational displacement (euler angles)
INTEGER(IntKi) :: I !< Generic counter

REAL(ReKi) :: outputArray(13)

CHARACTER(47) :: PointsOutputFmt !< Format specifier for the output file for wave elevation series
CHARACTER(3) :: AngleUnit !< Units for the angle

Expand Down Expand Up @@ -1972,32 +1973,30 @@ SUBROUTINE PointsForce_OutputWrite(ProgInfo, OutUnit, OutFileName, InputFileName
WRITE (OutUnit,'(A)', IOSTAT=ErrStatTmp ) ''
ENDIF


rotdisp = GetSmllRotAngs ( u%PtfmMesh%Orientation(:,:,1), ErrStatTmp, ErrMsgTmp )
CALL SetErrStat( ErrStatTmp, ErrMsgTmp, ErrStat, ErrMsg, RoutineName )


IF ( AnglesInDegrees ) THEN
CALL WrNumAryFileNR( OutUnit, (/ REAL(Time, ReKi), &
u%PtfmMesh%TranslationDisp(1,1), u%PtfmMesh%TranslationDisp(2,1), u%PtfmMesh%TranslationDisp(3,1), &
rotdisp(1)*R2D, rotdisp(2)*R2D, rotdisp(3)*R2D, &
u%PtfmMesh%TranslationVel(1,1), u%PtfmMesh%TranslationVel(2,1), u%PtfmMesh%TranslationVel(3,1), &
u%PtfmMesh%RotationVel(1,1)*R2D, u%PtfmMesh%RotationVel(2,1)*R2D, u%PtfmMesh%RotationVel(3,1)*R2D /), &
'3x,ES10.3E2', ErrStatTmp, ErrMsgTmp )
CALL WrNumAryFileNR( OutUnit, y%WriteOutput, '3x,ES10.3E2', ErrStatTmp, ErrMsgTmp )
WRITE (OutUnit,'(A)', IOSTAT=ErrStatTmp ) ''
outputArray = (/ REAL(Time, ReKi), &
REAL(u%PtfmMesh%TranslationDisp(1,1), ReKi), &
REAL(u%PtfmMesh%TranslationDisp(2,1), ReKi), &
REAL(u%PtfmMesh%TranslationDisp(3,1), ReKi), &
rotdisp(1)*R2D, rotdisp(2)*R2D, rotdisp(3)*R2D, &
u%PtfmMesh%TranslationVel(1,1), u%PtfmMesh%TranslationVel(2,1), u%PtfmMesh%TranslationVel(3,1), &
u%PtfmMesh%RotationVel(1,1)*R2D, u%PtfmMesh%RotationVel(2,1)*R2D, u%PtfmMesh%RotationVel(3,1)*R2D /)
ELSE
CALL WrNumAryFileNR( OutUnit, (/ REAL(Time, ReKi), &
u%PtfmMesh%TranslationDisp(1,1), u%PtfmMesh%TranslationDisp(2,1), u%PtfmMesh%TranslationDisp(3,1), &
rotdisp(1), rotdisp(2), rotdisp(3), &
u%PtfmMesh%TranslationVel(1,1), u%PtfmMesh%TranslationVel(2,1), u%PtfmMesh%TranslationVel(3,1), &
u%PtfmMesh%RotationVel(1,1), u%PtfmMesh%RotationVel(2,1), u%PtfmMesh%RotationVel(3,1) /), &
'3x,ES10.3E2', ErrStatTmp, ErrMsgTmp )
CALL WrNumAryFileNR( OutUnit, y%WriteOutput, '3x,ES10.3E2', ErrStatTmp, ErrMsgTmp )
WRITE (OutUnit,'(A)', IOSTAT=ErrStatTmp ) ''
outputArray = (/ REAL(Time, ReKi), &
REAL(u%PtfmMesh%TranslationDisp(1,1), ReKi), &
REAL(u%PtfmMesh%TranslationDisp(2,1), ReKi), &
REAL(u%PtfmMesh%TranslationDisp(3,1), ReKi), &
rotdisp(1), rotdisp(2), rotdisp(3), &
u%PtfmMesh%TranslationVel(1,1), u%PtfmMesh%TranslationVel(2,1), u%PtfmMesh%TranslationVel(3,1), &
u%PtfmMesh%RotationVel(1,1), u%PtfmMesh%RotationVel(2,1), u%PtfmMesh%RotationVel(3,1) /)
ENDIF



CALL WrNumAryFileNR( OutUnit, outputArray, '3x,ES10.3E2', ErrStatTmp, ErrMsgTmp )
CALL WrNumAryFileNR( OutUnit, y%WriteOutput, '3x,ES10.3E2', ErrStatTmp, ErrMsgTmp )
WRITE (OutUnit,'(A)', IOSTAT=ErrStatTmp ) ''

END SUBROUTINE PointsForce_OutputWrite

Expand Down