Skip to content

Regarding the issue of stiffness matrix coordinate transformation: #64

@TianyuanWangi

Description

@TianyuanWangi

Description

Describe the bug here

Hi,

I am trying to obtain the mooring line stiffness matrix using MoorPy and transform it to a coupled point. The MoorDyn file is for the 15 MW wind turbine (https://github.com/IEAWindSystems/IEA-15-240-RWT/blob/master/OpenFAST/IEA-15-240-RWT-UMaineSemi/IEA-15-240-RWT-UMaineSemi_MoorDyn.dat), where the coupling reference point is the PRP, located at (0,0,0) under hydrostatic conditions. However, when I want to transform the reference point to the COG (0,0,-14.4), I am unsure whether it is possible to convert the mooring stiffness matrix to the COG using a transformation method?

Method 1:

rRPRToCOG                    =  [0;0;-14.4];
stiffPRPToCOGMoorpy          =  OWStransMatRefPoint([0;0;0], rRPRToCOG, stiffPRPMoorpy(1:6,1:6));

stiffPRPToCOGMoorpy =

   1.0e+08 *

    0.0007         0   -0.0000         0    0.0218         0
         0    0.0007         0   -0.0218         0    0.0000
   -0.0000         0    0.0006         0   -0.0000         0
         0   -0.0218         0    3.0647         0    0.0000
    0.0218         0   -0.0000         0    3.0647         0
         0    0.0000         0   -0.0000         0    2.5229
function [ skewMatrix ] =  OWScalcSkewMatrix( refPoint, targetPoint )
% Produces a skew matrix. Multiplying this array by a vector is equivalent to the cross produce of vector and that vector
% posVec is a vector from the reference point to the target point

posVec                  =  targetPoint - refPoint;

skewMatrix              =  [ 0         -posVec(3)  posVec(2);
                             posVec(3)  0         -posVec(1);
                            -posVec(2)  posVec(1)  0  ];

end
function [ targetMatrix ] =  OWStransMatRefPoint( refPoint, targetPoint, refMatrix )

skewMatrix                 =  OWScalcSkewMatrix(refPoint, targetPoint); % refPoint and targetPoint are reshaped in OWScalcSkewMatrix
transformMatrix            =  [ eye(3)    skewMatrix;
                                zeros(3)  eye(3)];

targetMatrix           =  transformMatrix'*refMatrix*transformMatrix;

end

%% ================================================================
Method 2
Another approach I considered is to obtain the stiffness at each fairlead and then transform it to COG. However, the results always differ from those obtained through the coordinate transformation method, and I am not sure where the issue lies.

initDisp                     =  [0 0 0 0 0 0]';
fileName                     =  './PRP/lines.txt';
[stiffPRPMoorpy, posPRPMoorpy, loadPRPMoorpy] =  OWSstiffMatrixMoorpy(fileName, initDisp);

rPRPToFairlead               =  [-58 0 -14; 29 50.229 -14; 29 -50.229 -14]';
rRPRToCOG                    =  [0;0;-14.4];

stiffFairToPRPMoorpy         =  zeros(6);
stiffFairToCOGMoorpy         =  zeros(6);
stiffFairToTBMoorpy          =  zeros(6);

for iLine                     =  1:size(rPRPToFairlead,2)
    index                     =  iLine*6 + [1:6];
    iStiffFairMoorpy          =  stiffPRPMoorpy(index,1:6);

    index                     =  6 + (iLine-1)*3 + [1:3];
    Sf                        =  OWScalcSkewMatrix([0;0;0], loadPRPMoorpy(index,:));
    Sr                        =  OWScalcSkewMatrix([0;0;0], -rPRPToFairlead(:,iLine));
    iStiffFairMoorpy(4:6,4:6) =  Sf*Sr;
    stiffFairToPRPMoorpy      =  stiffFairToPRPMoorpy + OWStransMatRefPoint([0;0;0], -rPRPToFairlead(:,iLine), iStiffFairMoorpy);

    Sr                        =  OWScalcSkewMatrix([0;0;0], -rPRPToFairlead(:,iLine)+rRPRToCOG);
    iStiffFairMoorpy(4:6,4:6) =  Sf*Sr;
    stiffFairToCOGMoorpy      =  stiffFairToCOGMoorpy + OWStransMatRefPoint([0;0;0], -rPRPToFairlead(:,iLine)+rRPRToCOG, iStiffFairMoorpy);

end

stiffFairToCOGMoorpy =

   1.0e+08 *

    0.0007         0   -0.0000         0    0.0218         0
         0    0.0007         0   -0.0218         0    0.0000
   -0.0000         0    0.0006         0   -0.0000         0
         0   -0.0218         0    2.1888         0    0.0000
    0.0218         0   -0.0000         0    2.1888         0
         0    0.0000         0   -0.0000         0    2.5229
function [ stiffness, position, load, tension ] =  OWSstiffMatrixMoorpy( fileName, initDisp )

if nargin                    <  2
    initDisp                 =  [0;0;0;0;0;0];
end

ms                           =  mp.System(pyargs('file', fileName));
ms.initialize();
body                         =  ms.bodyList{1};
body.setPosition(initDisp);
ms.solveEquilibrium();

stiffnessPython              =  ms.getSystemStiffnessA(pyargs('DOFtype', 'coupled', 'lines_only', 'true', 'rho', waterDensity, 'g', gravity));
stiffness                    =  double(py.array.array('d', py.numpy.nditer(stiffnessPython)));
stiffness                    =  reshape(stiffness, stiffnessPython.shape{1}, stiffnessPython.shape{2})';

positionPython               =  ms.getPositions(pyargs('DOFtype', 'coupled'));
position                     =  double(py.array.array('d', py.numpy.nditer(positionPython)));
position                     =  reshape(position, [length(position) 1]);   

loadPython                   =  ms.getForces(pyargs('DOFtype', 'coupled', 'lines_only', 'true'));
load                         =  double(py.array.array('d', py.numpy.nditer(loadPython)));
load                         =  reshape(load, [length(load) 1]);     

tensionPython                =  ms.getTensions();
tension                      =  double(py.array.array('d', py.numpy.nditer(tensionPython)));
tension                      =  reshape(tension, [length(tension) 1]);    

counter                      =  0;
for iPoint                   =  1:length(ms.pointList)
    point                    =  ms.pointList{iPoint};
    if double(point.attachedEndB) ==  1  
        counter              =  counter + 1;
        fairStiff            =  zeros(6);
        iStiffnessPython     =  point.getStiffnessA();
        iStiffness           =  double(py.array.array('d', py.numpy.nditer(iStiffnessPython)));
        iStiffness           =  reshape(iStiffness, iStiffnessPython.shape{1}, iStiffnessPython.shape{2})';
        fairStiff(1:3,1:3)   =  iStiffness;
        stiffness            =  [stiffness; fairStiff];

        iPositionGlo         =  double(py.array.array('d', point.r));
        position             =  [position; reshape(iPositionGlo, [length(iPositionGlo) 1])];
        iPositionLoc         =  double(py.array.array('d', body.rPointRel{counter}));
        position             =  [position; reshape(iPositionLoc, [length(iPositionLoc) 1])];

        iLoadPython          =  point.getForces();   
        iLoad                =  double(py.array.array('d', py.numpy.nditer(iLoadPython)));
        load                 =  [load; reshape(iLoad, [length(iLoad) 1])];
    end
end
end

%% =============================
Method 2.1
The result below is the same as the transformation method

initDisp                     =  [0 0 0 0 0 0]';
fileName                     =  './PRP/lines.txt';
[stiffPRPMoorpy, posPRPMoorpy, loadPRPMoorpy] =  OWSstiffMatrixMoorpy(fileName, initDisp);

rPRPToFairlead               =  [-58 0 -14; 29 50.229 -14; 29 -50.229 -14]';
rRPRToCOG                    =  [0;0;-14.4];

stiffFairToPRPMoorpy         =  zeros(6);
stiffFairToCOGMoorpy         =  zeros(6);
stiffFairToTBMoorpy          =  zeros(6);

for iLine                     =  1:size(rPRPToFairlead,2)
    index                     =  iLine*6 + [1:6];
    iStiffFairMoorpy          =  stiffPRPMoorpy(index,1:6);

    index                     =  6 + (iLine-1)*3 + [1:3];
    Sf                        =  OWScalcSkewMatrix([0;0;0], loadPRPMoorpy(index,:));
    Sr                        =  OWScalcSkewMatrix([0;0;0], -rPRPToFairlead(:,iLine));
    iStiffFairMoorpy(4:6,4:6) =  Sf*Sr;
    stiffFairToPRPMoorpy      =  stiffFairToPRPMoorpy + OWStransMatRefPoint([0;0;0], -rPRPToFairlead(:,iLine), iStiffFairMoorpy);

    % Sr                        =  OWScalcSkewMatrix([0;0;0], -rPRPToFairlead(:,iLine)+rRPRToCOG);
    % iStiffFairMoorpy(4:6,4:6) =  Sf*Sr;
    stiffFairToCOGMoorpy      =  stiffFairToCOGMoorpy + OWStransMatRefPoint([0;0;0], -rPRPToFairlead(:,iLine)+rRPRToCOG, iStiffFairMoorpy);

end

stiffFairToCOGMoorpy =

   1.0e+08 *

    0.0007         0   -0.0000         0    0.0218         0
         0    0.0007         0   -0.0218         0    0.0000
   -0.0000         0    0.0006         0   -0.0000         0
         0   -0.0218         0    3.0647         0    0.0000
    0.0218         0   -0.0000         0    3.0647         0
         0    0.0000         0   -0.0000         0    2.5229

%% ===============================================
Method 3:

I also have an alternative modeling approach: defining the MoorDyn relative reference point as a vector relative to the COG and then obtaining the stiffness at the COG by specifying an initial displacement. This modeling approach produces results consistent with the second method.

initDisp                    =  [0 0 -14.4 0 0 0]';
fileName                    =  './COGInitDisp/lines.txt';
[stiffCOGInitDispMoorpy, posCOGInitDispMoorpy, loadCOGInitDispMoorpy] =  OWSstiffMatrixMoorpy(fileName, initDisp);

----------------------------------------------------------- MoorDyn INPUT FILE
IEA 15 MW offshore reference model on UChaine VolturnUS-S semi-submersible floating platform
----------------------------------------------------------- LINE TYPES
Name           Diam           MassDen        EA             BA/-zeta       EI             Cd             Ca             CdAx           CaAx
(-)            (m)            (kg/m)         (N)            (N-s/-)        (-)            (-)            (-)            (-)            (-)
main           0.333          685.000        3.270E+09     -1.000          0.000          2.000          0.820          0.400          0.270
----------------------------------------------------------- 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              0              0           0              0              0              0.0            0              0               0               0              0
----------------------------------------------------------- POINTS                                           
ID             Attachment     X              Y              Z              M              V              CdA            CA
(-)            (-)            (m)            (m)            (m)            (kg)           (m^3)          (m^2)          (-)
1              Body1         -58.000         0.000          0.400          0.000          0.000          0.000          0.000
2              Fixed         -837.600        0.000         -200.000        0.000          0.000          0.000          0.000
3              Body1          29.000         50.229         0.400          0.000          0.000          0.000          0.000
4              Fixed          418.800        725.383       -200.000        0.000          0.000          0.000          0.000
5              Body1          29.000        -50.229         0.400          0.000          0.000          0.000          0.000
6              Fixed          418.800       -725.383       -200.000        0.000          0.000          0.000          0.000
----------------------------------------------------------- LINES                                         
ID             LineType       AttachA        AttachB        UnstrLen       NumSegs        Outputs
(-)            (-)       	  (-)            (-)            (m)            (-)            (-)
1              main           2              1              850.000        50             -
2              main           4              3              850.000        50             -
3              main           6              5              850.000        50             -
----------------------------------------------------------- SOLVER OPTIONS
0.001                                                       dtM             - time step to use in mooring integration (s)
3.000E+06                                                   kbot            - bottom stiffness (Pa/m)
3.000E+05                                                   cbot      		- bottom damping (Pa-s/m)
1.000                                                       dtIC      		- time interval for analyzing convergence during IC gen (s)
120.000                                                     TmaxIC    		- max time for ic gen (s)
4.000                                                       CdScaleIC 		- factor by which to scale drag coefficients during dynamic relaxation (-)
0.001                                                       threshIC  		- threshold for IC convergence (-)    
RK4                                                         tScheme         - time integrator (RK4 is available for moorDyn v2)
200                                                         WtrDpth         - water depth
9.80665                                                     g               - gravity acceleration (m/s^2)
1025                                                        rho             - water density (kg/m^3)
0                                                           WaveKin         - wave kinematics flag (0=neglect, the only option currently supported)
0                                                           writeLog        - write a log file
0                                                           WriteUnits      - write the units header on the output files (must be zero, or report an error)  
1                                                           disableOutput   - disables some console and file outputs to improve runtime
1                                                           disableOutTime  - disables the printing of the current timestep to the console, useful for running with MATLAB                   
----------------------------------------------------------- OUTPUTS                                  
END of input file (the word "END" must appear in the first 3 columns of this last OutList line)
-----------------------------------------------------------                                             

Thank you very much for your advice

Steps to reproduce issue

Please provide a minimum working example (MWE) if possible

Current behavior

Expected behavior

Code versions

List versions only if relevant

  • Python

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions