Merge pull request #3321 from andrew-platt/b/SlD-SD_react_loads_backport3319

Backport #3319 SoilDyn/SubDyn: reaction loads when REDWIN is used

Author: andrew-platt

modules/openfast-library/src/FAST_Subs.f90

@@ -907,6 +907,9 @@ SUBROUTINE FAST_InitializeAll( t_initial, m_Glue, p_FAST, y_FAST, m_FAST, ED, SE
                     SlD%y, SlD%m, dt_module, Init%OutData_SlD, ErrStat2, ErrMsg2)
       if (Failed()) return
 
+      ! Pass nonlinear flag to SubDyn: true only when REDWIN DLL is active (CalcOption=3)
+      Init%InData_SD%SlDNonLinear = SlD%p%UseREDWINinterface ! REDWIN DLL returning nonlinear soil reaction forces      
+
       ! Add module to list of modules, return on error
       CALL MV_AddModule(m_Glue%ModData, Module_SlD, 'SlD', 1, dt_module, p_FAST%DT, &
                         Init%OutData_SlD%Vars, p_FAST%Linearize, ErrStat2, ErrMsg2)

modules/soildyn/src/SoilDyn.f90

@@ -109,9 +109,11 @@ subroutine SlD_Init(InitInp, u, p, x, xd, z, OtherState, y, m, Interval, InitOut
    z%DummyConstrState = 0.0_ReKi
    OtherState%DummyOtherState = 0.0_ReKi
 
-   ! are the returned reaction forces only the non-linear portion (used when SubDyn is calculating the linear portion)
+   ! The linear soil stiffness is passed directly to SubDyn as a stiffness matrix (not as a load)
+   ! The nonlinear portion of reaction forces is only computed when REDWIN DLL is active and passed as loads.
    p%SlDNonLinearForcePortionOnly = InitInp%SlDNonLinearForcePortionOnly
-   if (p%SlDNonLinearForcePortionOnly) call WrScr(' SoilDyn returning only non-linear portion of reaction forces')
+   call WrScr(' SoilDyn: linear stiffness passed directly to SubDyn as a stiffness matrix')
+   if (p%CalcOption == Calc_REDWIN) call WrScr(' SoilDyn: nonlinear soil reaction loads computed by REDWIN')
 
    ! If the module does not implement the four Jacobian routines at the end of this template, or the module cannot
    ! linearize with the features that are enabled, stop the simulation if InitInp%Linearize is true.

modules/subdyn/src/SubDyn.f90

@@ -479,6 +479,9 @@ SUBROUTINE SD_Init( InitInput, u, p, x, xd, z, OtherState, y, m, Interval, InitO
    ! Initialize the outputs & Store mapping between nodes and elements  
    CALL SDOUT_Init( Init, y, p, m, InitOut, InitInput%WtrDpth, ErrStat2, ErrMsg2 ); if(Failed()) return
 
+   ! Flag from glue code: if SoilDyn is returning nonlinear loads
+   p%SlDNonLinear = InitInput%SlDNonLinear
+
    ! Determine if we need to perform output file handling
    IF ( p%OutSwtch == 1 .OR. p%OutSwtch == 3 ) THEN  
        CALL SDOUT_OpenOutput( SD_ProgDesc, Init%RootName, p, InitOut, ErrStat2, ErrMsg2 ); if(Failed()) return

modules/subdyn/src/SubDyn_Output.f90

@@ -443,39 +443,45 @@ SUBROUTINE SDOut_MapOutputs(u,p,x, y, m, AllOuts, ErrStat, ErrMsg )
    ! --------------------------------------------------------------------------------{
    ! Total base reaction forces and moments at the (0.,0.,-WtrDpth) location in SS coordinate system
    !    "ReactFXss, ReactFYss, ReactFZss, ReactMXss, ReactMYss, ReactMZss"
-   IF (p%OutReact) THEN 
-      ALLOCATE ( ReactNs(6*p%nNodes_C), STAT = ErrStat )
-      IF ( ErrStat /= ErrID_None ) THEN
-         ErrMsg  = ' Error allocating space for ReactNs array.'
-         ErrStat = ErrID_Fatal
-         RETURN
-      END IF
-      ReactNs = 0.0_ReKi !Initialize
-      DO I=1,p%nNodes_C   !Do for each constrained node, they are ordered as given in the input file and so as in the order of y2mesh
-         FK_elm2=0._ReKi !Initialize for cumulative force
-         FM_elm2=0._ReKi !Initialize
-         pLst => p%MOutLst3(I)
-         !Find the joint forces
-         DO J=1,SIZE(pLst%ElmIDs(1,:))  !for all the elements connected (normally 1)
-            iiNode = 1
-            call ElementForce(pLst, iiNode, J, FM_elm, FK_elm, sgn, DIRCOS, .false.)
-            !transform back to global, need to do 3 at a time since cosine matrix is 3x3
-            DO L=1,2  
-               FM_elm2((L-1)*3+1:L*3) = FM_elm2((L-1)*3+1:L*3) + matmul(transpose(DIRCOS),FM_elm((L-1)*3+1:L*3))  !sum forces at joint in GLOBAL REF
-               FK_elm2((L-1)*3+1:L*3) = FK_elm2((L-1)*3+1:L*3) + matmul(transpose(DIRCOS),FK_elm((L-1)*3+1:L*3))  !signs may be wrong, we will fix that later;  
-               ! I believe this is all fixed in terms of signs now ,RRD 5/20/13
-            ENDDO           
+   IF (p%OutReact) THEN
+      IF (p%SlDNonLinear) THEN
+         ! When SoilDyn nonlinear loads are active (e.g., SoilDyn CalcOption = 3), SubDyn reaction loads are incomplete (only the linear part is included)
+         ! The total reaction at each base reaction joint is available in the SoilDyn output sensors (e.g., "Sld1Fxg Sld1Fyg Sld1Fzg Sld1Mxg Sld1Myg Sld1Mzg")
+         AllOuts( ReactSS(1:nDOFL_TP) ) = NaN
+      ELSE
+         ALLOCATE ( ReactNs(6*p%nNodes_C), STAT = ErrStat )
+         IF ( ErrStat /= ErrID_None ) THEN
+            ErrMsg  = ' Error allocating space for ReactNs array.'
+            ErrStat = ErrID_Fatal
+            RETURN
+         END IF
+         ReactNs = 0.0_ReKi !Initialize
+         DO I=1,p%nNodes_C   !Do for each constrained node, they are ordered as given in the input file and so as in the order of y2mesh
+            FK_elm2=0._ReKi !Initialize for cumulative force
+            FM_elm2=0._ReKi !Initialize
+            pLst => p%MOutLst3(I)
+            !Find the joint forces
+            DO J=1,SIZE(pLst%ElmIDs(1,:))  !for all the elements connected (normally 1)
+               iiNode = 1
+               call ElementForce(pLst, iiNode, J, FM_elm, FK_elm, sgn, DIRCOS, .false.)
+               !transform back to global, need to do 3 at a time since cosine matrix is 3x3
+               DO L=1,2
+                  FM_elm2((L-1)*3+1:L*3) = FM_elm2((L-1)*3+1:L*3) + matmul(transpose(DIRCOS),FM_elm((L-1)*3+1:L*3))  !sum forces at joint in GLOBAL REF
+                  FK_elm2((L-1)*3+1:L*3) = FK_elm2((L-1)*3+1:L*3) + matmul(transpose(DIRCOS),FK_elm((L-1)*3+1:L*3))  !signs may be wrong, we will fix that later;
+                  ! I believe this is all fixed in terms of signs now ,RRD 5/20/13
+               ENDDO
+            ENDDO
+            ! FK_elm2 ! + FM_elm2  !removed the inertial component 12/13 !Not sure why I need an intermediate step here, but the sum would not work otherwise
+            ! NEED TO ADD HYDRODYNAMIC FORCES AT THE RESTRAINT NODES
+            iSDNode   = p%Nodes_C(I,1)
+            iMeshNode = iSDNode ! input and Y2 mesh nodes are the same as subdyn
+            Fext =  (/ u%LMesh%Force(:,iMeshNode), u%LMesh%Moment(:,iMeshNode) /) + p%FG(p%NodesDOF(iMeshNode)%List(1:6))
+            Fext(1:3) = Fext(1:3) + p%FC(p%NodesDOF(iMeshNode)%List(1:3))
+            ReactNs((I-1)*6+1:6*I) = FK_elm2 - Fext  !Accumulate reactions from all nodes in GLOBAL COORDINATES
          ENDDO
-         ! FK_elm2 ! + FM_elm2  !removed the inertial component 12/13 !Not sure why I need an intermediate step here, but the sum would not work otherwise
-         ! NEED TO ADD HYDRODYNAMIC FORCES AT THE RESTRAINT NODES
-         iSDNode   = p%Nodes_C(I,1) 
-         iMeshNode = iSDNode ! input and Y2 mesh nodes are the same as subdyn
-         Fext =  (/ u%LMesh%Force(:,iMeshNode), u%LMesh%Moment(:,iMeshNode) /) + p%FG(p%NodesDOF(iMeshNode)%List(1:6))
-         Fext(1:3) = Fext(1:3) + p%FC(p%NodesDOF(iMeshNode)%List(1:3))
-         ReactNs((I-1)*6+1:6*I) = FK_elm2 - Fext  !Accumulate reactions from all nodes in GLOBAL COORDINATES
-      ENDDO
-      ! Store into AllOuts
-      AllOuts( ReactSS(1:nDOFL_TP) ) = matmul(p%TIreact,ReactNs)
+         ! Store into AllOuts
+         AllOuts( ReactSS(1:nDOFL_TP) ) = matmul(p%TIreact,ReactNs)
+      ENDIF
    ENDIF
    if (allocated(ReactNs)) deallocate(ReactNs)
 contains

modules/subdyn/src/SubDyn_Registry.txt

@@ -78,6 +78,7 @@ typedef  ^  InitInputType  ReKi             TP_RefPoint  {3}{:} - -  "global pos
 typedef  ^  InitInputType  ReKi             SubRotateZ   - - -       "Rotation angle in degrees about global Z"
 typedef  ^  InitInputType  ReKi             SoilStiffness ::: - -    "Soil stiffness matrices from SoilDyn" '(N/m, N-m/rad)'
 typedef  ^  InitInputType  MeshType         SoilMesh     - - -       "Mesh for soil stiffness locations" -
+typedef  ^  InitInputType  Logical          SlDNonLinear - - -       "Flag indicating that SoilDyn is returning nonlinear loads"
 typedef  ^  InitInputType  Logical          Linearize    - .FALSE. - "Flag that tells this module if the glue code wants to linearize." -
 
 # ============================== Initialization outputs ============================================================================================================================================
@@ -208,6 +209,7 @@ typedef  ^  ParameterType  IntKi   nDOFM        -  -  -  "retained degrees of fr
 typedef  ^  ParameterType  IntKi   nDOFRB       -  -  -  "number of rigid-body modes"
 typedef  ^  ParameterType  IntKi   SttcSolve    -  -  -  "Solve dynamics about static equilibrium point (flag)"
 typedef  ^  ParameterType  Logical Floating     -  -  -  "True if floating bottom (the 6 DOF are free at all reaction nodes)"
+typedef  ^  ParameterType  Logical SlDNonLinear -  -  -  "Flag indicating that SoilDyn is returning nonlinear loads"
 typedef  ^  ParameterType  ReKi    KMMDiag  {:}    -  -  "Diagonal coefficients of Kmm (OmegaM squared)"
 typedef  ^  ParameterType  ReKi    CMMDiag  {:}    -  -  "Diagonal coefficients of Cmm (~2 Zeta OmegaM))"
 typedef  ^  ParameterType  ReKi    MMB  {:}{:}     -  -  "Matrix after C-B reduction (transpose of MBM"

modules/subdyn/src/SubDyn_Types.f90

@@ -115,6 +115,7 @@ MODULE SubDyn_Types
     REAL(ReKi)  :: SubRotateZ = 0.0_ReKi      !< Rotation angle in degrees about global Z [-]
     REAL(ReKi) , DIMENSION(:,:,:), ALLOCATABLE  :: SoilStiffness      !< Soil stiffness matrices from SoilDyn ['(N/m,]
     TYPE(MeshType)  :: SoilMesh      !< Mesh for soil stiffness locations [-]
+    LOGICAL  :: SlDNonLinear = .false.      !< Flag indicating that SoilDyn is returning nonlinear loads [-]
     LOGICAL  :: Linearize = .FALSE.      !< Flag that tells this module if the glue code wants to linearize. [-]
   END TYPE SD_InitInputType
 ! =======================
@@ -258,6 +259,7 @@ MODULE SubDyn_Types
     INTEGER(IntKi)  :: nDOFRB = 0_IntKi      !< number of rigid-body modes [-]
     INTEGER(IntKi)  :: SttcSolve = 0_IntKi      !< Solve dynamics about static equilibrium point (flag) [-]
     LOGICAL  :: Floating = .false.      !< True if floating bottom (the 6 DOF are free at all reaction nodes) [-]
+    LOGICAL  :: SlDNonLinear = .false.      !< Flag indicating that SoilDyn is returning nonlinear loads [-]
     REAL(ReKi) , DIMENSION(:), ALLOCATABLE  :: KMMDiag      !< Diagonal coefficients of Kmm (OmegaM squared) [-]
     REAL(ReKi) , DIMENSION(:), ALLOCATABLE  :: CMMDiag      !< Diagonal coefficients of Cmm (~2 Zeta OmegaM)) [-]
     REAL(ReKi) , DIMENSION(:,:), ALLOCATABLE  :: MMB      !< Matrix after C-B reduction (transpose of MBM [-]
@@ -987,6 +989,7 @@ subroutine SD_CopyInitInput(SrcInitInputData, DstInitInputData, CtrlCode, ErrSta
    call MeshCopy(SrcInitInputData%SoilMesh, DstInitInputData%SoilMesh, CtrlCode, ErrStat2, ErrMsg2 )
    call SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
    if (ErrStat >= AbortErrLev) return
+   DstInitInputData%SlDNonLinear = SrcInitInputData%SlDNonLinear
    DstInitInputData%Linearize = SrcInitInputData%Linearize
 end subroutine
 
@@ -1023,6 +1026,7 @@ subroutine SD_PackInitInput(RF, Indata)
    call RegPack(RF, InData%SubRotateZ)
    call RegPackAlloc(RF, InData%SoilStiffness)
    call MeshPack(RF, InData%SoilMesh) 
+   call RegPack(RF, InData%SlDNonLinear)
    call RegPack(RF, InData%Linearize)
    if (RegCheckErr(RF, RoutineName)) return
 end subroutine
@@ -1044,6 +1048,7 @@ subroutine SD_UnPackInitInput(RF, OutData)
    call RegUnpack(RF, OutData%SubRotateZ); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpackAlloc(RF, OutData%SoilStiffness); if (RegCheckErr(RF, RoutineName)) return
    call MeshUnpack(RF, OutData%SoilMesh) ! SoilMesh 
+   call RegUnpack(RF, OutData%SlDNonLinear); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpack(RF, OutData%Linearize); if (RegCheckErr(RF, RoutineName)) return
 end subroutine
 
@@ -2413,6 +2418,7 @@ subroutine SD_CopyParam(SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg)
    DstParamData%nDOFRB = SrcParamData%nDOFRB
    DstParamData%SttcSolve = SrcParamData%SttcSolve
    DstParamData%Floating = SrcParamData%Floating
+   DstParamData%SlDNonLinear = SrcParamData%SlDNonLinear
    if (allocated(SrcParamData%KMMDiag)) then
       LB(1:1) = lbound(SrcParamData%KMMDiag)
       UB(1:1) = ubound(SrcParamData%KMMDiag)
@@ -3466,6 +3472,7 @@ subroutine SD_PackParam(RF, Indata)
    call RegPack(RF, InData%nDOFRB)
    call RegPack(RF, InData%SttcSolve)
    call RegPack(RF, InData%Floating)
+   call RegPack(RF, InData%SlDNonLinear)
    call RegPackAlloc(RF, InData%KMMDiag)
    call RegPackAlloc(RF, InData%CMMDiag)
    call RegPackAlloc(RF, InData%MMB)
@@ -3666,6 +3673,7 @@ subroutine SD_UnPackParam(RF, OutData)
    call RegUnpack(RF, OutData%nDOFRB); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpack(RF, OutData%SttcSolve); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpack(RF, OutData%Floating); if (RegCheckErr(RF, RoutineName)) return
+   call RegUnpack(RF, OutData%SlDNonLinear); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpackAlloc(RF, OutData%KMMDiag); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpackAlloc(RF, OutData%CMMDiag); if (RegCheckErr(RF, RoutineName)) return
    call RegUnpackAlloc(RF, OutData%MMB); if (RegCheckErr(RF, RoutineName)) return