search

NREL / MoorPy

BSD 3-Clause "New" or "Revised" License
28 stars 15 forks source link

Change in stiffness matrix when selecting plots parameter in solveEquilibrium #3

Closed Ahirvoas closed 2 years ago

Ahirvoas commented 2 years ago

Dear Moorpy developers,

I noticed a strange behavior of MoorPy 0.9.0 when estimating the stiffness matrix of the OC3-Hywind system. Indeed, the value of the lastly mentioned matrix varies depending on the selection (or not) of the plots parameter in the solveEquilibrium method.

Below is the input file I used with MoorPy:

--------------------- MoorDyn Input File ------------------------------------ Mooring system for OC3-Hywind TRUE Echo - echo the input file data (flag) ---------------------- LINE DICTIONARY ----------------------------------------------------- LineType Diam MassDenInAir EA cIntDamp EI Can Cat Cdn Cdt (-) (m) (kg/m) (N) (Pa-s) (N-m^2) (-) (-) (-) (-) main 0.09 77.7066 384.243E6 -0.8 1.0 0.0 1.6 0.1 --------------------- ROD DICTIONARY ----------------------------------------------------- RodType Diam MassDenInAir Can Cat Cdn Cdt (-) (m) (kg/m) (-) (-) (-) (-)
----------------------- BODY LIST ----------------------------------- BodyID X0 Y0 Z0 r0 p0 y0 Xcg Ycg Zcg M V IX IY IZ CdA-x,y,z Ca-x,y,z (-) (m) (m) (m) (deg) (deg) (deg) (m) (m) (m) (kg) (m^3) (kg-m^2) (kg-m^2) (kg-m^2) (m^2) (-) 1 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 -89.9155 7466.33E3 8026.5 4229.23E6 4229.23E6 164.23E6 0 0 0 0 ---------------------- ROD LIST -------------------- RodID Type/BodyID RodType Xa Ya Za Xb Yb Zb NumSegs Flags/Outputs (-) (-) (-) (m) (m) (m) (m) (m) (m) (-) (-)
---------------------- POINT LIST ----------------------------------------------------- Node Type X Y Z M V FX FY FZ CdA Ca (-) (-) (m) (m) (m) (kg) (m^3) (kN) (kN) (kN) (m2) () 1 Vessel 5.2 0.0 -70.0 0 0 0 0 0 0 0 2 Vessel -2.6 4.5 -70.0 0 0 0 0 0 0 0 3 Vessel -2.6 -4.5 -70.0 0 0 0 0 0 0 0 4 Fixed 853.87 0.0 -320.0 0 0 0 0 0 0 0 5 Fixed -426.94 739.47 -320.0 0 0 0 0 0 0 0 6 Fixed -426.94 -739.47 -320.0 0 0 0 0 0 0 0 ---------------------- LINE LIST ----------------------------------------------------- Line LineType UnstrLen NumSegs NodeAnch NodeFair Flags/Outputs (-) (-) (m) (-) (-) (-) (-) 1 main 902.2 20 4 1 pt 2 main 902.2 20 5 2 pt 3 main 902.2 20 6 3 pt ---------------------- SOLVER OPTIONS --------------------------------------- 0.001 dtM - time step to use in mooring integration 0 WaveKin - wave kinematics flag (1=include(unsupported), 0=neglect, 3=currentprofile.txt) 3.0e+06 kb - bottom stiffness 3.0e+05 cb - bottom damping 320.00 WtrDpth - water depth 4.0 ICDfac - factor by which to scale drag coefficients during dynamic relaxation IC gen 0.001 ICthresh - threshold for IC convergence 60 ICTmax - threshold for IC convergence ------------------------ OUTPUTS -------------------------------------------- FairTen1 FairTen2 FairTen3 AnchTen1 AnchTen2 AnchTen3 END ------------------------- need this line --------------------------------------

And hereafter, the python script used to compute the stiffness matrices:

import numpy as np
import moorpy as mp
# Load the MorrDyn file
ms = mp.System('NRELOffshrBsline5MW_OC3Hywind_MoorDyn.dat', depth=320)

ms.initialize(plots=0)
ms.solveEquilibrium(DOFtype="free", plots=0, rmsTol=10, maxIter=200)
K_mooring  = ms.getSystemStiffness(DOFtype="free", dx = 0.1, dth = 0.1, solveOption=1, plots=0)
ms_1 = mp.System('NRELOffshrBsline5MW_OC3Hywind_MoorDyn.dat', depth=320)
ms_1.initialize(plots=0)
ms_1.solveEquilibrium(DOFtype="free", plots=1, rmsTol=10, maxIter=200)
K_mooring_1  = ms_1.getSystemStiffness(DOFtype="free", dx = 0.1, dth = 0.1, solveOption=1, plots=0)

Then, I figured out, thanks to the source code, that the method getSystemStiffness is calling solveEquilibrium function with plots parameter sets to zero. Thus, I am now using the following approach to get the stiffness matrix:

import numpy as np
import moorpy as mp
ms = mp.System('NRELOffshrBsline5MW_OC3Hywind_MoorDyn.dat', depth=320)
ms.initialize(plots=0)
K_mooring = ms.getSystemStiffness(DOFtype="free", dx = 0.1, dth = 0.1, solveOption=1, plots=0)

Best regards,

shousner commented 2 years ago

Hi Ahirvoas,

Thanks for reaching out. It's very exciting to see the code being used by more people! Apologies for the delayed response.

In general, MoorPy is still very much in a "development" phase and the documentation is still very much in progress.

It seems like your approach and logic are all correct, but I'll offer some suggestions to get things working:

ms = mp.System('NRELOffshrBsline5MW_OC3Hywind_MoorDyn.dat', depth=320)
ms.initialize(plots=0)
print(ms.bodyList[0].type)

If the above printed line is 0, then the body is a free body. If it's 1, it's fixed. If it's -1, it's coupled

ms.solveEquilibrium(DOFtype="free", plots=0, rmsTol=10, maxIter=200)

This call to solveEquilibrium (see the recent commit) will not change anything if the body type is 1 or -1, since there are no free points to equilibrate If the body type is 0, then the body will want to balance forces, which include weight/buoyancy (see the Body class)

Now that the system is equilibrated, we can calculate the stiffness matrix If the body type is 0, I would recommend

K_mooring1 = ms.getSystemStiffness(DOFtype="free", dx = 0.1, dth = 0.1, solveOption=1, plots=0)

If the body is -1, I would recommend

K_mooring2 = ms.getCoupledStiffness(any optional inputs you want)

Regardless of the body type, I would still recommend running

K_mooringA = ms.getSystemStiffnessA(options)

And comparing this to any of the other system stiffness matrices, especially in terms of run time (it's a much more efficient method)

Hopefully this clears up any issues between solveEquilibrium and any system stiffness method.

Stein

mattEhall commented 2 years ago

Presumably resolve, so closing this issue.