Dear Jason,
Thanks for your prompt reply. Thanks too for your suggestion. The code I snipped and inserted into my post above was from TwrLoading. I passed TwrAM and TwrFt to the new subroutine (which calculates external loads on the tower due to aerodynamic drag) because that was how MorisonTwrLd and UserTwrLd were set up. I will change the code so that a local array is used to store an intemediate value of the external tower load and then copied to FTAero (containing only the aerodynamic force per unit length).
Kind Regards,
Mark
[code]SUBROUTINE TwrLoading ( JNode, X1 , X2 , X3 , X4 , X5 , X6 , &
XD1, XD2, XD3, XD4, XD5, XD6 )
! This routine computes the tower hydrodynamic loading; that is
! TwrAM(1:6,1:6) and TwrFt(1:6).
USE General
USE FixedBottomSupportStructure, ONLY:MorisonTwrLd
USE SimCont
USE Tower
USE Features, ONLY: CompAero
USE RtHndSid, ONLY: FTAero, MFAero
IMPLICIT NONE
! Passed Variables:
REAL(ReKi), INTENT(IN ) :: X1 ! The xi-component of the translational displacement (in m ) of the current tower node relative to the inertial frame origin at ground level [onshore] or MSL [offshore].
REAL(ReKi), INTENT(IN ) :: X2 ! The yi-component of the translational displacement (in m ) of the current tower node relative to the inertial frame origin at ground level [onshore] or MSL [offshore].
REAL(ReKi), INTENT(IN ) :: X3 ! The zi-component of the translational displacement (in m ) of the current tower node relative to the inertial frame origin at ground level [onshore] or MSL [offshore].
REAL(ReKi), INTENT(IN ) :: X4 ! The xi-component of the rotational displacement (in rad ) of the current tower element relative to the inertial frame origin at ground level [onshore] or MSL [offshore].
REAL(ReKi), INTENT(IN ) :: X5 ! The yi-component of the rotational displacement (in rad ) of the current tower element relative to the inertial frame origin at ground level [onshore] or MSL [offshore].
REAL(ReKi), INTENT(IN ) :: X6 ! The zi-component of the rotational displacement (in rad ) of the current tower element relative to the inertial frame origin at ground level [onshore] or MSL [offshore].
REAL(ReKi), INTENT(IN ) :: XD1 ! The xi-component of the translational velocity (in m/s ) of the current tower node relative to the inertial frame origin at ground level [onshore] or MSL [offshore].
REAL(ReKi), INTENT(IN ) :: XD2 ! The yi-component of the translational velocity (in m/s ) of the current tower node relative to the inertial frame origin at ground level [onshore] or MSL [offshore].
REAL(ReKi), INTENT(IN ) :: XD3 ! The zi-component of the translational velocity (in m/s ) of the current tower node relative to the inertial frame origin at ground level [onshore] or MSL [offshore].
REAL(ReKi), INTENT(IN ) :: XD4 ! The xi-component of the rotational (angular) velocity (in rad/s) of the current tower element relative to the inertial frame origin at ground level [onshore] or MSL [offshore].
REAL(ReKi), INTENT(IN ) :: XD5 ! The yi-component of the rotational (angular) velocity (in rad/s) of the current tower element relative to the inertial frame origin at ground level [onshore] or MSL [offshore].
REAL(ReKi), INTENT(IN ) :: XD6 ! The zi-component of the rotational (angular) velocity (in rad/s) of the current tower element relative to the inertial frame origin at ground level [onshore] or MSL [offshore].
INTEGER(4), INTENT(IN ) :: JNode ! The number of the current tower node / element. [1 to TwrNodes]
! Local variables:
REAL(ReKi), PARAMETER :: SymTol = 9.999E-4 ! Tolerance used to determine if matrix PtfmAM is symmetric.
REAL(ReKi) :: X (6) ! The 3 components of the translational displacement (in m ) of the current tower node and the 3 components of the rotational displacement (in rad ) of the current tower element relative to the inertial frame origin at ground level [onshore] or MSL [offshore].
REAL(ReKi) :: XD (6) ! The 3 components of the translational velocity (in m/s) of the current tower node and the 3 components of the rotational (angular) velocity (in rad/s) of the current tower element relative to the inertial frame origin at ground level [onshore] or MSL [offshore].
REAL(ReKi) :: LocTwrFrc(6) ! local copy of tower forces to ensure that hydrodynamic effects are superposed on top of aerodynamic effects
INTEGER(4) :: I ! Loops through all platform DOFs.
INTEGER(4) :: J ! Loops through all platform DOFs.
! Place the displacement and velocity arguments into the local arrays,
! X(1:6) and XD(1:6), respectively:
X (1) = X1
X (2) = X2
X (3) = X3
X (4) = X4
X (5) = X5
X (6) = X6
XD(1) = XD1
XD(2) = XD2
XD(3) = XD3
XD(4) = XD4
XD(5) = XD5
XD(6) = XD6
! Compute the tower aerodynamic loading for the current tower node/ element:
IF ( CompAero ) THEN ! AeroDyn has been used
CALL CalcAeroDrag ( JNode, X, XD, ZTime, DirRoot, TwrAM, LocTwrFrc )
FTAero(JNode,1) = LocTwrFrc(1) !
FTAero(JNode,2) = LocTwrFrc(2) !
FTAero(JNode,3) = LocTwrFrc(3) ! write output tower forces and moments into
MFAero(JNode,1) = LocTwrFrc(4) ! FTAero and MFAero
MFAero(JNode,2) = LocTwrFrc(5) !
MFAero(JNode,3) = LocTwrFrc(6) !
ENDIF
! Compute the tower hydrodynamic loading for the current tower node/ element:
SELECT CASE ( PtfmModel ) ! Which platform model are we using?
CASE ( 0 ) ! None!
! user has selected not to calculate loads on the platform.
! no calculation required here, since TwrAM and TwrFt are all initialized to zero.
CASE ( 1 ) ! Onshore.
! user has selected not to calculate loads on the platform.
! no calculation required here, since TwrAM and TwrFt are all initialized to zero.
CASE ( 2 ) ! Fixed bottom offshore.
SELECT CASE ( TwrLdMod ) ! Which tower loading model are we using?
CASE ( 0 ) ! None!
! user has selected not to calculate loads on the tower.
! no calculation required here, since TwrAM and TwrFt are all initialized to zero.
CASE ( 1 ) ! Undocumented hydrodynamic loading using Morison’s equation.
! CALL the undocumented Morison's equation tower loading model:
LocTwrFrc(1) = TwrFt(1)
LocTwrFrc(2) = TwrFt(2)
LocTwrFrc(3) = TwrFt(3)
CALL MorisonTwrLd ( JNode, DiamT(JNode), CAT(JNode), CDT(JNode), X, XD, ZTime, TwrAM, TwrFt )
IF ( ( TwrFt(1) /= 0 .AND. LocTwrFrc(1) /= 0 ) .OR. ( TwrFt(2) /= 0 .AND. LocTwrFrc(2) /= 0 ) .OR. ( TwrFt(3) /= 0 .AND. LocTwrFrc(3) /= 0 ) ) THEN
CALL ProgWarn( ' Warning: both aerodynamic tower loads (CalcAeroDrag) and hydrodynamic loads (MorisonTwrLd) are simultaneously non-zero.' )
ENDIF
TwrFt(1) = TwrFt(1) + LocTwrFrc(1)
TwrFt(2) = TwrFt(2) + LocTwrFrc(2)
TwrFt(3) = TwrFt(3) + LocTwrFrc(3)
CASE ( 2 ) ! User-defined tower loading.
! CALL the user-defined tower loading model:
CALL UserTwrLd ( JNode, X, XD, ZTime, DirRoot, TwrAM, TwrFt )
! Ensure that the tower element added mass matrix returned by UserTwrLd,
! TwrAM, is symmetric; Abort if necessary:
DO I = 1,5 ! Loop through the 1st 5 rows (columns) of TwrAM
DO J = (I+1),6 ! Loop through all columns (rows) passed I
IF ( ABS( TwrAM(I,J) - TwrAM(J,I) ) > SymTol ) &
CALL ProgAbort ( ' The user-defined tower element added mass matrix is unsymmetric.'// &
' Make sure TwrAM returned by UserTwrLd() is symmetric.' )
ENDDO ! J - All columns (rows) passed I
ENDDO ! I - The 1st 5 rows (columns) of TwrAM
ENDSELECT
CASE ( 3 ) ! Floating offshore.
! Do nothing here since TwrAM and TwrFt are all initialized to zero.
ENDSELECT
RETURN
END SUBROUTINE TwrLoading
!=======================================================================
SUBROUTINE CalcAeroDrag ( JNode, X, XD, ZTime, DirRoot, TwrAM, TwrFt )
! CalcAeroDrag has been written to calculate aerodynamic drag loads
! on the tower (the portion of the support structure above the current
! water level and below the yaw bearing).
!
! The tower loads returned by this routine contain contributions from
! any external load acting on the current tower element (indicated by
! JNode) other than loads transmitted from tower aerodynamics.
!
! No contribution to the local tower mass has been calculated below.
!
! This routine was originally written assuming that the tower loads are
! transmitted through a medium like soil [foundation] and/or water
! [offshore], so that added mass effects are important. Consequently,
! the routine assumes that the total load per unit length on the current
! tower element can be written as:
!
! TwrF(i) = SUM( -TwrAM(i,j)*XDD(j), j=1,2,…,6) + TwrFt(i) for i=1,2,…,6
!
! where,
! TwrF(i) = the i’th component of the total load per unit length
! applied on the current tower element; positive in the
! direction of positive motion of the i’th DOF of the current
! tower element
! TwrAM(i,j) = the (i,j) component of the tower added mass matrix per unit
! length (output by this routine)
! XDD(j) = the j’th component of the current tower element
! acceleration vector
! TwrFt(i) = the i’th component of the portion of the current tower
! element load per unit length associated with everything but
! the added mass effects; positive in the direction of
! positive motion of the i’th DOF of the current tower
! element (output by this routine)
! The order of indices in all arrays passed to and from this routine is as
! follows:
! 1 = Current tower element surge / xi-component of translation
! 3 = Current tower element sway / yi-component of translation
! 3 = Current tower element heave / zi-component of translation
! 4 = Current tower element roll / xi-component of rotation
! 5 = Current tower element pitch / yi-component of rotation
! 6 = Current tower element yaw / zi-component of rotation
! NOTE: The added mass matrix returned by this routine, TwrAM, must be
! symmetric. FAST and ADAMS will abort otherwise.
!
! Please also note that the hydrostatic restoring contribution to the
! hydrodynamic force returned by this routine should not contain the
! effects of body weight, as is often done in classical marine
! hydrodynamics. The effects of body weight are included within FAST
! and ADAMS.
USE AeroDyn, ONLY: TwrAeroDrag
USE Tower, ONLY: TwrNodes
USE TurbConf, ONLY: TowerHt, TwrRBHT, TwrDraft
USE Waves, ONLY: WaveElevation, WaveKinzi0, DZNodes, WaveStMod
IMPLICIT NONE
! Passed Variables:
REAL(ReKi), INTENT(OUT) :: TwrAM (6,6) ! Added mass matrix per unit length of current tower element, kg/m, kg-m/m, kg-m^2/m.
REAL(ReKi), INTENT(OUT) :: TwrFt (6) ! The surge/xi (1), sway/yi (2), and heave/zi (3)-components of the portion of the tower force per unit length (in N/m) at the current tower element and the roll/xi (4), pitch/yi (5), and yaw/zi (6)-components of the portion of the tower moment per unit length (in N-m/m) acting at the current tower element associated with everything but the added-mass effects; positive forces are in the direction of motion.
REAL(ReKi), INTENT(IN ) :: X (6) ! The 3 components of the translational displacement (in m ) of the current tower node and the 3 components of the rotational displacement (in rad ) of the current tower element relative to the inertial frame origin at ground level [onshore] or MSL [offshore].
REAL(ReKi), INTENT(IN ) :: XD (6) ! The 3 components of the translational velocity (in m/s) of the current tower node and the 3 components of the rotational (angular) velocity (in rad/s) of the current tower element relative to the inertial frame origin at ground level [onshore] or MSL [offshore].
REAL(ReKi), INTENT(IN ) :: ZTime ! Current simulation time, sec.
INTEGER, INTENT(IN ) :: JNode ! The number of the current tower node / element, (-). [1 to TwrNodes]
CHARACTER(1024), INTENT(IN ) :: DirRoot ! The name of the root file including the full path to the current working directory. This may be useful if you want this routine to write a permanent record of what it does to be stored with the simulation results: the results should be stored in a file whose name (including path) is generated by appending any suitable extension to DirRoot.
! Local Variables:
REAL(ReKi) :: Damp (6,6) ! Damping matrix.
REAL(ReKi) :: Stff (6,6) ! Stiffness/restoring matrix.
REAL(ReKi), DIMENSION(1:3) :: TwrVel ! The surge/xi (1), sway/yi (2), and heave/zi (3)-components of local tower linear velocity [m/s]
REAL(ReKi), DIMENSION(1:3) :: TwrFrc ! The surge/xi (1), sway/yi (2), and heave/zi (3)-components of the external drag force on the tower node, JNode
REAL(ReKi), DIMENSION(1:3) :: TwrPos ! The surge/xi (1), sway/yi (2), and heave/zi (3)-components of the current deflection of the tower node, JNode
REAL(ReKi), DIMENSION(1:3) :: TwrGeom ! The three components of tower geometry required by TwrAeroDrag (TwrRBHt, TwrDraft and TowerHt)
INTEGER(4) :: I ! Generic index.
INTEGER(4) :: J ! Generic index.
LOGICAL :: OK2CalcAeroDrag ! .TRUE. if entire element is above the water level (including the height of any incident wave)
REAL(ReKi) :: WaveElevation0 ! current local value of wave elevation
TwrPos(1) = X(1)
TwrPos(2) = X(2)
TwrPos(3) = X(3)
TwrFrc = 0.0
! calculate aerodynamic drag on tower if the local tower element mid-point is above the instantaneous water height
!
! some definitions:-
!
! (1) WaveStMod – Model for stretching incident wave kinematics to instantaneous free surface
! {0: none=no stretching, 1: vertical stretching, 2: extrapolation stretching, 3: Wheeler stretching}
! (2) WaveKinzi0(1:NWaveKin0) – zi-coordinates for points along a vertical line passing through the platform reference point
! where the incident wave kinematics will be computed; these are relative to the mean see
! level (meters)
! (3) DZNodes(1:NWaveKin0) – Length of variable-length tower or platform elements for the points along a vertical line
! passing through the platform reference point where the incident wave kinematics will be
! computed (meters)
! (4) WaveElevation0 – Elevation of incident waves at the platform reference point and current time (meters)
! .TRUE. if we have no stretching; therefore, integrate up to the MSL, regardless of the instantaneous free surface elevation.
IF ( WaveStMod == 0 ) THEN
! .TRUE. if the current tower element lies entirely above the MSL.
IF ( ( ( WaveKinzi0(JNode) - 0.5DZNodes(JNode) ) >= 0.0 ) .AND. ( TwrPos(3) - 0.5DZNodes(JNode) > 0.0 ) ) THEN
OK2CalcAeroDrag = .TRUE.
ELSE
OK2CalcAeroDrag = .FALSE.
ENDIF
ELSE ! We must have some sort of stretching.
WaveElevation0 = WaveElevation ( 1, ZTime ) ! compute the local wave elevation
! .TRUE. if the current tower element lies entirely above the free surface of the incident wave.
IF ( ( ( WaveKinzi0(JNode) - 0.5DZNodes(JNode) ) >= WaveElevation0 ) .AND. ( TwrPos(3) - 0.5DZNodes(JNode) > 0.0 ) ) THEN
OK2CalcAeroDrag = .TRUE.
ELSE
OK2CalcAeroDrag = .FALSE.
ENDIF
ENDIF
IF ( OK2CalcAeroDrag ) THEN
TwrGeom(1) = TwrRBHt
TwrGeom(2) = TwrDraft
TwrGeom(3) = TowerHt
TwrVel(1) = XD(1)
TwrVel(2) = XD(2)
TwrVel(3) = XD(3)
! get tower local diameter and drag coefficient
CALL TwrAeroDrag(ZTime, TwrPos, TwrVel, TwrFrc, JNode, TwrNodes, TwrGeom)
ELSE
TwrFrc(1) = 0.0
TwrFrc(2) = 0.0
TwrFrc(3) = 0.0
ENDIF
Damp (1,
= (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
Damp (2, = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
Damp (3, = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
Damp (4, = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
Damp (5, = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
Damp (6, = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
Stff (1, = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
Stff (2, = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
Stff (3, = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
Stff (4, = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
Stff (5, = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
Stff (6, = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
TwrAM(1, = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
TwrAM(2, = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
TwrAM(3, = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
TwrAM(4, = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
TwrAM(5, = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
TwrAM(6, = (/ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 /)
TwrFt(1) = TwrFrc(1)
TwrFt(2) = TwrFrc(2)
TwrFt(3) = TwrFrc(3)
TwrFt(4) = 0.0 ! no bending moments applied due to aerodynamic drag on the tower
TwrFt(5) = 0.0
TwrFt(6) = 0.0
DO J = 1,6
DO I = 1,6
TwrFt(I) = TwrFt(I) - Damp(I,J)*XD(J) - Stff(I,J)*X(J)
ENDDO
ENDDO
RETURN
END SUBROUTINE CalcAeroDrag
!=======================================================================
END MODULE FASTSubs
[/code]
