Skip to content

MD: Adding Load dependent dynamic stiffness #2459

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

Merged
andrew-platt merged 12 commits into from
Oct 21, 2024
Merged

MD: Adding Load dependent dynamic stiffness #2459

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
71c0a2f
Updated r-tests
RyanDavies19 Aug 9, 2024
d8ca05b
Update r-test pointer
RyanDavies19 Aug 9, 2024
f4ba9ee
Update r-test pointer
RyanDavies19 Sep 27, 2024
76dad55
Merge remote-tracking branch 'upstream/dev' into dev
RyanDavies19 Sep 27, 2024
4cbcb9a
Updated copyright, adjusted driver statement, updated dynamic stiffness
RyanDavies19 Oct 2, 2024
5ab5eaa
Viscoelastic bug fixes: state allocation and assinging variables
RyanDavies19 Oct 11, 2024
92d3c8f
Update viscoelastic varible names, fix dynamic stiffness model
RyanDavies19 Oct 11, 2024
7743f6f
Add viscoelastic error checking
RyanDavies19 Oct 11, 2024
2fb06e3
Merge remote-tracking branch 'upstream/dev' into update_EA_D
RyanDavies19 Oct 11, 2024
36020b7
typo fix in MD viscoelastic
RyanDavies19 Oct 17, 2024
bab4d45
One more typo fix
RyanDavies19 Oct 17, 2024
a15e465
Revert incidental changes to Registr_NWTC_Library
andrew-platt Oct 21, 2024
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
23 changes: 15 additions & 8 deletions modules/moordyn/src/MoorDyn.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -167,7 +167,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er
InitOut%Ver = MD_ProgDesc

CALL WrScr(' This is MoorDyn v2, with significant input file changes from v1.')
CALL DispCopyrightLicense( MD_ProgDesc%Name, 'Copyright (C) 2019 Matt Hall' )
CALL DispCopyrightLicense( MD_ProgDesc%Name)


!---------------------------------------------------------------------------------------------
Expand Down Expand Up @@ -638,10 +638,14 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er

! process stiffness coefficients
CALL SplitByBars(tempString1, N, tempStrings)
if (N > 2) then
CALL SetErrStat( ErrID_Fatal, 'A line type EA entry can have at most 2 (comma-separated) values.', ErrStat, ErrMsg, RoutineName )
if (N > 3) then
CALL SetErrStat( ErrID_Fatal, 'A line type EA entry can have at most 3 (bar-separated) values.', ErrStat, ErrMsg, RoutineName )
CALL CleanUp()
else if (N==2) then ! visco-elastic case!
else if (N==3) then ! visco-elastic case, load dependent dynamic stiffness!
m%LineTypeList(l)%ElasticMod = 3
read(tempStrings(2), *) m%LineTypeList(l)%alphaMBL
read(tempStrings(3), *) m%LineTypeList(l)%vbeta
else if (N==2) then ! visco-elastic case, constant dynamic stiffness!
m%LineTypeList(l)%ElasticMod = 2
read(tempStrings(2), *) m%LineTypeList(l)%EA_D
else
Expand All @@ -657,11 +661,11 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er
! process damping coefficients
CALL SplitByBars(tempString2, N, tempStrings)
if (N > m%LineTypeList(l)%ElasticMod) then
CALL SetErrStat( ErrID_Fatal, 'A line type BA entry cannot have more (comma-separated) values its EA entry.', ErrStat, ErrMsg, RoutineName )
CALL SetErrStat( ErrID_Fatal, 'A line type BA entry cannot have more (bar-separated) values than its EA entry.', ErrStat, ErrMsg, RoutineName )
CALL CleanUp()
else if (N==2) then ! visco-elastic case when two BA values provided
read(tempStrings(2), *) m%LineTypeList(l)%BA_D
else if (m%LineTypeList(l)%ElasticMod == 2) then ! case where there is no dynamic damping for viscoelastic model (will it work)?
else if (m%LineTypeList(l)%ElasticMod > 1) then ! case where there is no dynamic damping for viscoelastic model (will it work)?
CALL WrScr("Warning, viscoelastic model being used with zero damping on the dynamic stiffness.")
if (p%writeLog > 0) then
write(p%UnLog,'(A)') "Warning, viscoelastic model being used with zero damping on the dynamic stiffness."
Expand Down Expand Up @@ -1438,7 +1442,7 @@ SUBROUTINE MD_Init(InitInp, u, p, x, xd, z, other, y, m, DTcoupling, InitOut, Er

! account for states of line
m%LineStateIs1(l) = Nx + 1
if (m%LineTypeList(m%LineList(l)%PropsIdNum)%ElasticMod == 2) then
if (m%LineTypeList(m%LineList(l)%PropsIdNum)%ElasticMod > 1) then ! todo add an error check here? or change to 2 or 3?
Nx = Nx + 7*m%LineList(l)%N - 6 ! if using viscoelastic model, need one more state per segment
m%LineStateIsN(l) = Nx
else
Expand Down Expand Up @@ -3498,7 +3502,10 @@ SUBROUTINE MD_CalcContStateDeriv( t, u, p, x, xd, z, other, m, dxdt, ErrStat, Er

! calculate line dynamics (and calculate line forces and masses attributed to points)
DO l = 1,p%nLines
CALL Line_GetStateDeriv(m%LineList(l), dxdt%states(m%LineStateIs1(l):m%LineStateIsN(l)), m, p) !dt might also be passed for fancy friction models
CALL Line_GetStateDeriv(m%LineList(l), dxdt%states(m%LineStateIs1(l):m%LineStateIsN(l)), m, p, ErrStat, ErrMsg) !dt might also be passed for fancy friction models
if (ErrStat == ErrID_Fatal) then
return
endif
END DO

! calculate point dynamics (including contributions from attached lines
Expand Down
2 changes: 1 addition & 1 deletion modules/moordyn/src/MoorDyn_Driver.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -144,7 +144,7 @@ PROGRAM MoorDyn_Driver
CALL CPU_TIME ( ProgStrtCPU ) ! Initial time (this zeros the start time when used as a MATLAB function)


CALL WrScr('MD Driver updated '//TRIM( version%Date ))
CALL WrScr('MD Driver last updated '//TRIM( version%Date ))

! Parse the driver input file and run the simulation based on that file
CALL get_command_argument(1, drvrFilename)
Expand Down
72 changes: 56 additions & 16 deletions modules/moordyn/src/MoorDyn_Line.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,11 +69,13 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg)
Line%d = LineProp%d
Line%rho = LineProp%w/(Pi/4.0 * Line%d * Line%d)

Line%EA = LineProp%EA
Line%EA = LineProp%EA
! note: Line%BA is set later
Line%EA_D = LineProp%EA_D
Line%BA_D = LineProp%BA_D
Line%EI = LineProp%EI !<<< for bending stiffness
Line%EA_D = LineProp%EA_D
Line%alphaMBL = LineProp%alphaMBL
Line%vbeta = LineProp%vbeta
Line%BA_D = LineProp%BA_D
Line%EI = LineProp%EI !<<< for bending stiffness

Line%Can = LineProp%Can
Line%Cat = LineProp%Cat
Expand All @@ -82,6 +84,12 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg)

! copy over elasticity data
Line%ElasticMod = LineProp%ElasticMod

if (Line%ElasticMod > 3) then
ErrStat = ErrID_Fatal
ErrMsg = "Line ElasticMod > 3. This is not possible."
RETURN
endif

Line%nEApoints = LineProp%nEApoints
DO I = 1,Line%nEApoints
Expand Down Expand Up @@ -141,7 +149,7 @@ SUBROUTINE SetupLine (Line, LineProp, p, ErrStat, ErrMsg)
END IF

! if using viscoelastic model, allocate additional state quantities
if (Line%ElasticMod == 2) then
if (Line%ElasticMod > 1) then
ALLOCATE ( Line%dl_1(N), STAT = ErrStat )
IF ( ErrStat /= ErrID_None ) THEN
ErrMsg = ' Error allocating dl_1 array.'
Expand Down Expand Up @@ -991,7 +999,7 @@ SUBROUTINE Line_SetState(Line, X, t)
END DO

! if using viscoelastic model, also set the static stiffness stretch
if (Line%ElasticMod == 2) then
if (Line%ElasticMod > 1) then
do I=1,Line%N
Line%dl_1(I) = X( 6*Line%N-6 + I) ! these will be the last N entries in the state vector
end do
Expand All @@ -1001,12 +1009,15 @@ END SUBROUTINE Line_SetState
!--------------------------------------------------------------

!--------------------------------------------------------------
SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p) !, FairFtot, FairMtot, AnchFtot, AnchMtot)
SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p, ErrStat, ErrMsg) !, FairFtot, FairMtot, AnchFtot, AnchMtot)

TYPE(MD_Line), INTENT(INOUT) :: Line ! the current Line object
Real(DbKi), INTENT(INOUT) :: Xd(:) ! state derivative vector section for this line
TYPE(MD_MiscVarType), INTENT(INOUT) :: m ! passing along all mooring objects
TYPE(MD_ParameterType), INTENT(IN ) :: p ! Parameters
INTEGER(IntKi), INTENT( OUT) :: ErrStat ! Error status of the operation
CHARACTER(*), INTENT( OUT) :: ErrMsg ! Error message if ErrStat /= ErrID_None


! Real(DbKi), INTENT( IN ) :: X(:) ! state vector, provided
! Real(DbKi), INTENT( INOUT ) :: Xd(:) ! derivative of state vector, returned ! cahnged to INOUT
Expand Down Expand Up @@ -1044,7 +1055,8 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p) !, FairFtot, FairMtot, AnchFtot,
Real(DbKi) :: Yi ! used in interpolating from lookup table
Real(DbKi) :: dl ! stretch of a segment [m]
Real(DbKi) :: ld_1 ! rate of change of static stiffness portion of segment [m/s]
Real(DbKi) :: EA_1 ! stiffness of 'static stiffness' portion of segment, combines with dynamic stiffness to give static stiffnes [m/s]
Real(DbKi) :: EA_1 ! stiffness of 'slow' portion of segment, combines with dynamic stiffness to give static stiffnes [m/s]
Real(DbKi) :: EA_D ! stiffness of 'fast' portion of segment, combines with EA_1 stiffness to give static stiffnes [m/s]

REAL(DbKi) :: surface_height ! Average the surface heights at the two nodes
REAL(DbKi) :: firstNodeZ ! Difference of first node depth from surface height
Expand Down Expand Up @@ -1268,21 +1280,49 @@ SUBROUTINE Line_GetStateDeriv(Line, Xd, m, p) !, FairFtot, FairMtot, AnchFtot,
else
MagT = 0.0_DbKi ! cable can't "push"
end if

! line internal damping force based on line-specific BA value, including possibility of dynamic length changes in l and ld terms
MagTd = Line%BA* ( Line%lstrd(I) - Line%lstr(I)*Line%ld(I)/Line%l(I) )/Line%l(I)

! viscoelastic model
else if (Line%ElasticMod == 2) then
! viscoelastic model from https://asmedigitalcollection.asme.org/OMAE/proceedings/IOWTC2023/87578/V001T01A029/1195018
else if (Line%ElasticMod > 1) then

if (Line%ElasticMod == 3) then
if (Line%dl_1(I) >= 0.0) then
! Mean load dependent dynamic stiffness: from combining eqn. 2 and eqn. 10 from original MD viscoelastic paper, taking mean load = k1 delta_L1 / MBL, and solving for k_D using WolframAlpha with following conditions: k_D > k_s, (MBL,alpha,beta,unstrLen,delta_L1) > 0
EA_D = 0.5 * ((Line%alphaMBL) + (Line%vbeta*Line%dl_1(I)*(Line%EA / Line%l(I))) + Line%EA + sqrt((Line%alphaMBL * Line%alphaMBL) + (2*Line%alphaMBL*(Line%EA / Line%l(I)) * (Line%vbeta*Line%dl_1(I) - Line%l(I))) + ((Line%EA / Line%l(I))*(Line%EA / Line%l(I)) * (Line%vbeta*Line%dl_1(I) + Line%l(I))*(Line%vbeta*Line%dl_1(I) + Line%l(I)))))
else
EA_D = Line%alphaMBL ! mean load is considered to be 0 in this case. The second term in the above equation is not valid for delta_L1 < 0.
endif

else if (Line%ElasticMod == 2) then
! constant dynamic stiffness
EA_D = Line%EA_D
endif

if (EA_D == 0.0) then ! Make sure EA != EA_D or else nans, also make sure EA_D != 0 or else nans.
ErrStat = ErrID_Fatal
ErrMsg = "Viscoelastic model: Dynamic stiffness cannot equal zero"
return
else if (EA_D == Line%EA) then
ErrStat = ErrID_Fatal
ErrMsg = "Viscoelastic model: Dynamic stiffness cannot equal static stiffness"
return
endif

EA_1 = Line%EA_D*Line%EA/(Line%EA_D - Line%EA)! calculated EA_1 which is the stiffness in series with EA_D that will result in the desired static stiffness of EA_S
EA_1 = EA_D*Line%EA/(EA_D - Line%EA)! calculated EA_1 which is the stiffness in series with EA_D that will result in the desired static stiffness of EA_S.

dl = Line%lstr(I) - Line%l(I) ! delta l of this segment

ld_1 = (Line%EA_D*dl - (Line%EA_D + EA_1)*Line%dl_1(I) + Line%BA_D*Line%lstrd(I)) /( Line%BA_D + Line%BA) ! rate of change of static stiffness portion [m/s]

!MagT = (Line%EA*Line%dl_S(I) + Line%BA*ld_S)/ Line%l(I) ! compute tension based on static portion (dynamic portion would give same)
MagT = EA_1*Line%dl_1(I)/ Line%l(I)
MagTd = Line%BA*ld_1 / Line%l(I)
ld_1 = (EA_D*dl - (EA_D + EA_1)*Line%dl_1(I) + Line%BA_D*Line%lstrd(I)) /( Line%BA_D + Line%BA) ! rate of change of static stiffness portion [m/s]

if (dl >= 0.0) then ! if both spring 1 (the spring dashpot in parallel) and the whole segment are not in compression
MagT = EA_1*Line%dl_1(I) / Line%l(I) ! compute tension based on static portion (dynamic portion would give same). See eqn. 14 in paper
else
MagT = 0.0_DbKi ! cable can't "push"
endif

MagTd = Line%BA*ld_1 / Line%l(I) ! compute tension based on static portion (dynamic portion would give same). See eqn. 14 in paper

! update state derivative for static stiffness stretch (last N entries in the state vector)
Xd( 6*N-6 + I) = ld_1
Expand Down
Loading