!======================================================================= SUBROUTINE UserTwrLd ( JNode, X, XD, ZTime, DirRoot, TwrAM, TwrFt ) ! This routine implements the distributed springs foundation model ! for the NREL 5MW Offshore Baseline Wind Turbine. A direct CALL ! to routine MorisonTwrLd() is made in order to obtain the ! hydrodynamic loads on the monopile between the mudline and the ! free surface. The distributed springs at each tower node are ! determined by interpolating into the array of discrete springs ! specified by Patrik Passon of the University of Stuttgart in ! Germany as documented in: ! "OC3-Derivation and Description of the Soil-Pile-Interaction Models.pdf" ! "OC3-Soil-Pile_Interaction_Model.xls" ! 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 ! 2 = 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 ! This routine was written by J. Jonkman of NREL/NWTC on August 28, ! 2006 for use in the IEA Annex XXIII OC3 studies. USE FixedBottomSupportStructure USE InterpSubs USE Precision USE Waves, ONLY:WtrDpth, NWaveKin0, WaveKinzi0, DZNodes 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(4), 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) :: DZFract ! The fraction of the current tower element that is below the seabed (0.0 <= DZFract <= 1.0): 0.0 = the element is entirely above the seabed, 1.0 = element is entirely below the seabed (-) REAL(ReKi), ALLOCATABLE :: Spring (:) ! Discrete springs as specified by Patrik Passon (N/m) REAL(ReKi), ALLOCATABLE,SAVE :: Stff (:) ! The values of Spring(:) interpolated to each of the undeflected tower node elevations (as stored in the WaveKinzi0(:) array) (N/m) REAL(ReKi), PARAMETER :: TwrCA = 1.0 ! Normalized hydrodynamic added mass coefficient in Morison's equation (-) REAL(ReKi), PARAMETER :: TwrCD = 1.0 ! Normalized hydrodynamic viscous drag coefficient in Morison's equation (-) REAL(ReKi), PARAMETER :: TwrDiam = 6.0 ! Tower diameter in Morison's equation (meters) REAL(ReKi), ALLOCATABLE :: zi (:) ! Elevations (w.r.t. MSL) were the discrete springs are specified by Patrik Passon (meters) INTEGER(4) :: J ! Generic index INTEGER(4) :: K ! Generic index INTEGER(4) :: LastInd = 1 ! Index into the arrays saved from the last call as a starting point for this call INTEGER(4) :: NSpring ! Number of discrete springs specifed by Patrik Passon (-) INTEGER(4) :: Sttus ! Status returned by an attempted allocation. LOGICAL(1), SAVE :: FirstPas = .TRUE. ! When .TRUE., indicates we're on the first pass ! Perform initialization calculations on first pass only: IF ( FirstPas ) THEN ! Define the distributed springs using the properties specified by Patrik ! Passon: NSpring = 37 ALLOCATE ( zi (NSpring ) , STAT=Sttus ) IF ( Sttus /= 0 ) THEN CALL Abort(' Error allocating memory for the zi array.') ENDIF ALLOCATE ( Spring(NSpring ) , STAT=Sttus ) IF ( Sttus /= 0 ) THEN CALL Abort(' Error allocating memory for the Spring array.') ENDIF zi ( 1) = -56.0 zi ( 2) = -55.0 zi ( 3) = -54.0 zi ( 4) = -53.0 zi ( 5) = -52.0 zi ( 6) = -51.0 zi ( 7) = -50.0 zi ( 8) = -49.0 zi ( 9) = -48.0 zi (10) = -47.0 zi (11) = -46.0 zi (12) = -45.0 zi (13) = -44.0 zi (14) = -43.0 zi (15) = -42.0 zi (16) = -41.0 zi (17) = -40.0 zi (18) = -39.0 zi (19) = -38.0 zi (20) = -37.0 zi (21) = -36.0 zi (22) = -35.0 zi (23) = -34.0 zi (24) = -33.0 zi (25) = -32.0 zi (26) = -31.0 zi (27) = -30.0 zi (28) = -29.0 zi (29) = -28.0 zi (30) = -27.0 zi (31) = -26.0 zi (32) = -25.0 zi (33) = -24.0 zi (34) = -23.0 zi (35) = -22.0 zi (36) = -21.0 zi (37) = -20.0 Spring( 1) = 595318000.0 Spring( 2) = 1165208000.0 Spring( 3) = 1129400000.0 Spring( 4) = 1095553000.0 Spring( 5) = 1059931000.0 Spring( 6) = 1024493000.0 Spring( 7) = 989209000.0 Spring( 8) = 953643000.0 Spring( 9) = 918718000.0 Spring(10) = 883287000.0 Spring(11) = 847803000.0 Spring(12) = 812541000.0 Spring(13) = 777187000.0 Spring(14) = 741870000.0 Spring(15) = 706616000.0 Spring(16) = 671440000.0 Spring(17) = 636229000.0 Spring(18) = 600957000.0 Spring(19) = 565919000.0 Spring(20) = 530470000.0 Spring(21) = 495081000.0 Spring(22) = 459574000.0 Spring(23) = 385327000.0 Spring(24) = 305479000.0 Spring(25) = 280059000.0 Spring(26) = 254125000.0 Spring(27) = 227500000.0 Spring(28) = 200112000.0 Spring(29) = 171927000.0 Spring(30) = 143115000.0 Spring(31) = 114173000.0 Spring(32) = 80184000.0 Spring(33) = 52237000.0 Spring(34) = 35561000.0 Spring(35) = 20912000.0 Spring(36) = 9000000.0 Spring(37) = 1156000.0 ! Interpolate to find the foundation stiffness, Stff(:), at each of the ! undeflected tower node elevations (as stored in the WaveKinzi0(:) ! array): ALLOCATE ( Stff (NWaveKin0) , STAT=Sttus ) IF ( Sttus /= 0 ) THEN CALL Abort(' Error allocating memory for the Stff array.') ENDIF DO J = 1,NWaveKin0 ! Loop through the tower nodes / elements IF ( WaveKinzi0(J) > -WtrDpth ) THEN ! .TRUE. if the current node lies above the seabed. Stff(J) = 0.0 ELSE ! The current tower node must lie below the seabed. Stff(J) = InterpStp( WaveKinzi0(J), zi(:), Spring(:), LastInd, NSpring ) ENDIF ENDDO ! J - Tower nodes / elements FirstPas = .FALSE. ! Don't enter here again! ENDIF ! First, initialize the added mass matrix per unit length of the current ! tower element, TwrAM, and the portion of the current tower element load ! per unit length associated with everything but the added mass effects, ! TwrFt, using the hydrodynamic effects via a direct CALL to routine ! MorisonTwrLd(): CALL MorisonTwrLd ( JNode, TwrDiam, TwrCA, TwrCD, X, XD, ZTime, TwrAM, TwrFt ) ! Find the fraction of the current tower element that is below the seabed: IF ( ( WaveKinzi0(JNode) - 0.5*DZNodes(JNode) ) >= -WtrDpth ) THEN ! .TRUE. if the current tower element lies entirely above the seabed. DZFract = 0.0 ELSEIF ( ( WaveKinzi0(JNode) + 0.5*DZNodes(JNode) ) <= -WtrDpth ) THEN ! .TRUE. if the current tower element lies entirely below the seabed. DZFract = 1.0 ELSE ! The seabed must fall somewhere along the current tower element; thus, interpolate. DZFract = ( -WtrDpth - ( WaveKinzi0(JNode) - 0.5*DZNodes(JNode) ) )/DZNodes(JNode) ENDIF ! Compute the foundation loads using uncoupled linear springs for the ! portion of the current tower element that lies below the seabed: IF ( DZFract > 0.0 ) THEN ! .TRUE. if a portion of the current tower element lies below the seabed. DO K = 1,2 ! Loop through the xi- (1) and yi- (2) directions TwrFt(K) = TwrFt(K) - Stff(JNode)*X(K)*DZFract ENDDO ! K - The xi- (1) and yi- (2) directions ENDIF RETURN END SUBROUTINE UserTwrLd !=======================================================================