MoorPy Linearization for BModes

Ian.Ammerman · February 27, 2025, 8:56pm

Hi,

I would like to use MoorPy to get the mooring stiffness matrix for a floating platform. I have put the mooring system into MoorPy and it initializes well:

I then go to linearize the system using the following python code:

print('===== Beginning Linearization =====')

# Initialize stiffness matrix
K = np.zeros((6, 6), dtype=np.float64)

# Set displacements
DOFchange = np.array([6,6,6,0.2,0.2,0.2])

# Force threshold for linearization
Fmin = 1 # N or Nm

# Loop over DOF and get stiffness matrix
for DOF in range(0,6):

    print(f"----- Perturbing DOF {DOF+1} -----")

    # Positive displaced position
    newPos1 = np.array([0, 0, 0, 0, 0, 0], dtype=np.float64)
    newPos1[DOF] = DOFchange[DOF]
    # print(newPos1)

    # Negative displaced position
    newPos2 = np.array([0, 0, 0, 0, 0, 0], dtype=np.float64)
    newPos2[DOF] = -DOFchange[DOF]
    # print(newPos2)

    # Set positive position
    ms.bodyList[0].setPosition(newPos1)

    # Get forces
    ms.solveEquilibrium()
    fpos = ms.bodyList[0].getForces()
    fpos[np.absolute(fpos) < Fmin] = 0

    # Set negative position
    ms.bodyList[0].setPosition(newPos2)

    # Get forces
    ms.solveEquilibrium()
    fneg = ms.bodyList[0].getForces()
    fneg[np.absolute(fneg) < Fmin] = 0

    # Compute stiffness
    df = fpos - fneg
    dd = newPos2[DOF] - newPos1[DOF]
    print("dF: [" + " ".join(f"{x:.2e}" for x in df) + "]")
    Kcol = np.divide(df,dd)
    print("Kcol: [" + " ".join(f"{x:.2e}" for x in df) + "]")
    # print(Kcol)

    # Append column onto stiffness matrix
    K[:, DOF] = Kcol

It seems to do OK, but no matter what I do, I cannot get a symmetric mooring stiffness matrix. I am implementing a catenary mooring system for a spar platform from this paper with the following parameters:

The mooring stiffness matrix I get is this:

A couple of notes on my mind:

Some of these values seem quite small, and in the python code I include some logic to get rid of REALLY small values (less than 1N), but my intuition is that this should be symmetric for the mooring system. Is this possibly a limitation in MoorPy, or maybe there is a better approach I could take?
The coordinate system is rotated from the typical OpenFAST, no? So in MoorPy, +x is into the wind, +z is up? I am planning to use this stiffness matrix in BModes to get floating modes and frequencies of a turbine, but I suspect I need to negate these terms to match BModes and OpenFAST sign convention.

Best,
Ian

Roger.Bergua · February 27, 2025, 9:54pm

Hi @Ian.Ammerman,

Recently I did a similar exercise in MoorPy. See for reference: SubDyn vs MoorPy: Differences in the linearized stiffness matrix for horizontal pretensioned lines · Issue #2637 · OpenFAST/openfast · GitHub

To get the linearized stiffness matrix in MoorPy, I would recommend you to use:

KsystemA = ms.getSystemStiffnessA(DOFtype=‘coupled’)

This corresponds to the analytical calculation of the stiffness matrix.

As commented in the Github issue linked above, the stiffness matrix is not symmetric if there is an initial cable pretension and an assymmetry in the mooring lines. I believe this non-symmetric behavior should appear in the lower-right quadrant of the stiffness matrix.

Regarding the coordinate system, my understanding (based on the tests that I performed) is that the coordinate system is the same as OpenFAST. This is x pointing in the downstream direction and z pointing upwards.

I hope that helps!

Ian.Ammerman · February 27, 2025, 10:15pm

Hi @Roger.Bergua,

Thank you for the suggestion! I re-ran and compared to my result and they were quite similar. The internal function yielded a truly symmetric matrix as I had expected but it’s nice to know I was on the right path.

My mooring system does have pre-tension, but is symmetric so I would expect a reasonably symmetrix matrix. My result was always *nearly symmetric.

As for the coordinate system - yes, I believe you are right. I went back through how I generated the MoorDynV2 file from MoorPy (very handy!) and I had mixed up my headings by 180.

Roger.Bergua · February 28, 2025, 1:51am

I’m glad to hear that it worked for you, @Ian.Ammerman. Looking at your stiffness matrix, it seems that the coefficients k46 and k64 were very different. Could you confirm if you got a symmetric stiffness matrix?

Ian.Ammerman · February 28, 2025, 2:12am

Hi @Roger.Bergua ,

Here is the updated one that I computed. I set all elements with abs(K) < 1 to zero to make it easier to read.

Looking closely, you’re right that the K46/K64 is still not symmetric. I’m not sure what could do that - dynamically it should be the same I think.

Ian.Ammerman · February 28, 2025, 2:13am

For reference, here is the generated MoorDynV2 file:

MoorDyn v2 Input File 
Generated by MoorPy
---------------------- LINE TYPES --------------------------------------------------
TypeName      Diam     Mass/m     EA     BA/-zeta     EI        Cd      Ca      CdAx    CaAx
(name)        (m)      (kg/m)     (N)    (N-s/-)    (N-m^2)     (-)     (-)     (-)     (-)
chain1        0.0900   245.40  3.840e+08 -1.000e+00 0.000e+00   1.333   1.000   0.64    0.50   
--------------------- ROD TYPES -----------------------------------------------------
TypeName      Diam     Mass/m    Cd     Ca      CdEnd    CaEnd
(name)        (m)      (kg/m)    (-)    (-)     (-)      (-)
----------------------- BODIES ------------------------------------------------------
ID   Attachment    X0     Y0     Z0     r0      p0     y0     Mass          CG*          I*      Volume   CdA*   Ca*
(#)     (-)        (m)    (m)    (m)   (deg)   (deg)  (deg)   (kg)          (m)         (kg-m^2)  (m^3)   (m^2)  (-)
1     coupled     0.00   0.00   0.00   0.00   0.00   0.00   0.0000e+00  0.00|0.00|0.00 0.000e+00   0.00   0.00  0.00
---------------------- RODS ---------------------------------------------------------
ID   RodType  Attachment  Xa    Ya    Za    Xb    Yb    Zb   NumSegs  RodOutputs
(#)  (name)    (#/key)    (m)   (m)   (m)   (m)   (m)   (m)  (-)       (-)
---------------------- POINTS -------------------------------------------------------
ID  Attachment     X       Y       Z           Mass  Volume  CdA    Ca
(#)   (-)         (m)     (m)     (m)          (kg)  (m^3)  (m^2)   (-)
1    Fixed       855.20     0.00  -320.00      0.00   0.00   0.00   0.00
2    Body1         6.50     0.00   -70.00      0.00   0.00   0.00   0.00
3    Fixed      -427.60   740.62  -320.00      0.00   0.00   0.00   0.00
4    Body1        -3.25     5.63   -70.00      0.00   0.00   0.00   0.00
5    Fixed      -427.60  -740.62  -320.00      0.00   0.00   0.00   0.00
6    Body1        -3.25    -5.63   -70.00      0.00   0.00   0.00   0.00
---------------------- LINES --------------------------------------------------------
ID    LineType      AttachA  AttachB  UnstrLen  NumSegs  LineOutputs
(#)    (name)        (#)      (#)       (m)       (-)     (-)
1    chain1            1       2      902.200     40       p
2    chain1            3       4      902.200     40       p
3    chain1            5       6      902.200     40       p
---------------------- OPTIONS ------------------------------------------------------
0.001            dtM
3000000.0        kb
300000.0         cb
60               TmaxIC
0.001            dtm
60               tmaxic
9.81             g
320.0            depth
1025.0           rho
----------------------- OUTPUTS -----------------------------------------------------
FairTen1
FairTen2
FairTen3
END
--------------------- need this line ------------------------------------------------

Roger.Bergua · February 28, 2025, 3:40pm

Thanks for the feedback, @Ian.Ammerman.

The matrix that you attach has a small asymmetry in the positions k46 vs k46: 4.1692E3 vs 2.5581E4. However, this asymmetry is significantly different than the one that you showed in the first matrix: 766.4764 vs -11451470 (?)

From the last matrix provided by MoorPy, I can see that the stiffness in roll (k44) and pitch (k55) is slightly different. A similar behavior can be observed between surge (k11) and sway (k22) diections. This indicates that in reality, there is a [very] small asymmetry in your system.

It seems that this very small asymmetry comes from the points defined in the space. I took the line aligned with the x-axis and projected it 120 and 240 degrees, increasing the number of decimals positions (i.e., more precision).

By doing so, I get a stiffness matrix that is symmetric (as expected for this axisymmetric arrangement of mooring lines). Moreover, the k46 and k64 coefficients are in reality 0.

I hope this helps!

Ian.Ammerman · February 28, 2025, 4:30pm

Hi @Roger.Bergua ,

Very interesting - I just projected them myself to confirm your numbers and I agree. I get the same matrix after adapting the MoorDyn file.

The MoorDynV2 file I was using was generated by MoorPy like this:

# -------- Code to Define System Manually -------- #
# ------------------------------------------------ #
# ----- choose some system geometry parameters -----

depth     = 320                             # water depth [m]
angles    = np.radians([0, 120, 240])      # line headings list [rad]
rAnchor   = 855.2                            # anchor radius/spacing [m]
zFair     = -70                             # fairlead z elevation [m]
rFair     = 6.5                              # fairlead radius [m]
lineLength= 902.2                            # line unstretched length [m]
typeName  = "chain1"                        # identifier string for the line type


# ----- set up the mooring system and floating body -----

# Create new MoorPy System and set its depth
ms = mp.System(depth=depth)

# add a line type
ms.setLineType(dnommm=90, material='chain', name=typeName)  # this would be 120 mm chain

# Add a free, body at [0,0,0] to the system (including some properties to make it hydrostatically stiff)
ms.addBody(0, np.zeros(6), m=0, v=1e3, rM=100, AWP=1e3)
ms.bodyList[0].type = -1

# For each line heading, set the anchor point, the fairlead point, and the line itself
for i, angle in enumerate(angles):

    # create end Points for the line
    ms.addPoint(1, [rAnchor*np.cos(angle), rAnchor*np.sin(angle), -depth])   # create anchor point (type 0, fixed)
    ms.addPoint(1, [  rFair*np.cos(angle),   rFair*np.sin(angle),  zFair])   # create fairlead point (type 0, fixed)

    # attach the fairlead Point to the Body (so it's fixed to the Body rather than the ground)
    ms.bodyList[0].attachPoint(2*i+2, [rFair*np.cos(angle), rFair*np.sin(angle), zFair])

    # add a Line going between the anchor and fairlead Points
    ms.addLine(lineLength, typeName, pointA=2*i+1, pointB=2*i+2)
    
# Load mooring system from file
# ms = mp.System(file='IEA-10.0-198-RWT-MoorDyn-v2.txt')
ms.initialize()                                             # make sure everything's connected
print(ms.bodyList[0].r6)
ms.solveEquilibrium()                                       # equilibrate
fig, ax = ms.plot()                                         # plot the system in original configuration
ms.unload("IEA-10.0-198-RWT-MoorDyn-v2.txt")                         # export to MD input file

I further modified the line types values (Diam, Mass/m, EA, etc.) myself because I was unsure how to specify these in the MoorPy generation code directly (is that currently implemented?)

Given what we just saw in terms of effect of the significant figures, I modifed system.py lines which print the x,y,z coordinates to the file to automatically print more significant figures to the file:

L.append("{:<4d} {:9} {:8.6f} {:8.6f} {:8.6f} {:9.2f} {:6.2f} {:6.2f} {:6.2f}".format(
                              point.number,point_type, point_pos[0],point_pos[1],point_pos[2], point.m, point.v, point.CdA, point.Ca))

This makes me wonder if the same should be done for the other sections in the file, and if there is a minimum number of figures to ensure convergence. What are your thoughts on this?

Roger.Bergua · February 28, 2025, 4:54pm

Hi @Ian.Ammerman,

I don’t define the MoorPy input file by means of a code. I write it manually. For a very complex mooring system, I guess I would code it…

Personally, I would not be concerned about the decimal positions. I would use 2 decimal positions for the locations (i.e., precision in the order of centimeters).

I just wanted to illustrate and explain what was the reason for the stiffness matrix not being symmetric in this case where I was expecting a symmetric matrix given the symmetry in the mooring set up.

Topic		Replies	Views
About the linear stiffness matrix of the anchor chain Structural Analysis	5	493	October 4, 2023
Mooring stiffness matrix Structural Analysis	7	6744	April 2, 2018
Linearized mooring line stiffness matrix Computer-Aided Engineering Software Tools	8	738	August 16, 2024
Mooring modelization on openfast Systems Engineering Software Tools	2	1257	July 20, 2021
Moordyn convergence problem consultation Computer-Aided Engineering Software Tools	2	31	January 9, 2025

MoorPy Linearization for BModes

Related topics