[manoj_nav]

Hi All

I am trying to use interDyMfoam solver for doing dynamic trim and sinkage simulation of a ship. Following is an extract from the log.interDyMfoam file.

6-DoF rigid body motion

Centre of rotation: (2.59 0 0.296454)

Centre of mass: (2.59 0 0.296454)

Orientation: (0.99971 0 -0.02408 0 1 0 0.02408 0 0.99971)

Linear velocity: (0 0 -0.0128308)

Angular velocity: (0 0.00336935 0)



What is this orientation tensor? How do we get trim (rotation about y axis) and sinkage ( translation about z axis) information out of it?

Regards,

Manoj

[davibarreira]

When you run your simulation, save the log file. In the log file you will have the information that you want, so it's just a matter of writing a little script to extract it. Here is a little bash to get the center of mass.

#!/bin/bash
grep 'Centre of mass' $1 | cut -d '(' -f 2 | tr -d ")" > centreMass_Extract
grep -e "^Time = " $1 | cut -d " " -f 3 > times_Extract

paste times_Extract centreMass_Extract > cmMotion_parcial
sed -e 's/ [ ]*/\t/g' cmMotion_parcial > cmMotion
rm cmMotion_parcial
rm times_Extract
rm centreMass_Extract

[manoj_nav]

Originally Posted by davibarreira

When you run your simulation, save the log file. In the log file you will have the information that you want, so it's just a matter of writing a little script to extract it. Here is a little bash to get the center of mass.

#!/bin/bash
grep 'Centre of mass' $1 | cut -d '(' -f 2 | tr -d ")" > centreMass_Extract
grep -e "^Time = " $1 | cut -d " " -f 3 > times_Extract

paste times_Extract centreMass_Extract > cmMotion_parcial
sed -e 's/ [ ]*/\t/g' cmMotion_parcial > cmMotion
rm cmMotion_parcial
rm times_Extract
rm centreMass_Extract

Thank you for the reply. It works very well.

But I am facing 2 issues .

1) Because I am having 3 nOuterCorrectors loops inside PIMPLE, so it seems 6DOF solvers also runs the same number of time for every time step. So in the cmmotion file, center of mass data is 3 times the time data. I am attaching the files with the message. If you can help me to correct it. I am new to writing these scripts.

2) What is the orientation tensor? There is 9 values in it. What are these values?

Thanks again.

Regards,

Manoj

[manoj_nav]

Originally Posted by manoj_nav

1) Because I am having 3 nOuterCorrectors loops inside PIMPLE, so it seems 6DOF solvers also runs the same number of time for every time step. So in the cmmotion file, center of mass data is 3 times the time data. I am attaching the files with the message. If you can help me to correct it. I am new to writing these scripts.

Got this working by including,

awk 'NR == 1 || NR % 3 == 0' centreMass_Extract > centreMass_Extract1


But I still dont understand what are these 9 values of orientation tensor..

[manoj_nav]

Originally Posted by manoj_nav

But I still dont understand what are these 9 values of orientation tensor..

Can anyone confirm if the orientation tensor is give by below matrix, where the rotations are about the , and axes with angles , and ?

[kmooney]

Originally Posted by manoj_nav

Can anyone confirm if the orientation tensor is give by below matrix, where the rotations are about the , and axes with angles , and ?

Looking at this:
http://planning.cs.uiuc.edu/node102.html

and this:
http://foam.sourceforge.net/docs/cpp/a09729_source.html

It indeed appears that that is the correct summation of the Rx Ry and Rz tensors.

[manoj_nav]

Thanks a lot for the reply.

Regards,

Manoj

[haze_1986]

For interDyMFoam, if you are using sixDoFRigidBodyMotionI.H:

inline Foam::tensor Foam::sixDoFRigidBodyMotion::rotationTensorX
(
    scalar phi
) const
{
    return tensor
    (
        1, 0, 0,
        0, Foam::cos(phi), -Foam::sin(phi),
        0, Foam::sin(phi), Foam::cos(phi)
    );
}


inline Foam::tensor Foam::sixDoFRigidBodyMotion::rotationTensorY
(
    scalar phi
) const
{
    return tensor
    (
        Foam::cos(phi), 0, Foam::sin(phi),
        0, 1, 0,
        -Foam::sin(phi), 0, Foam::cos(phi)
    );
}


inline Foam::tensor Foam::sixDoFRigidBodyMotion::rotationTensorZ
(
    scalar phi
) const
{
    return tensor
    (
        Foam::cos(phi), -Foam::sin(phi), 0,
        Foam::sin(phi), Foam::cos(phi), 0,
        0, 0, 1
    );
}


inline Foam::Tuple2<Foam::tensor, Foam::vector>
Foam::sixDoFRigidBodyMotion::rotate
(
    const tensor& Q0,
    const vector& pi0,
    const scalar deltaT
) const
{
    Tuple2<tensor, vector> Qpi(Q0, pi0);
    tensor& Q = Qpi.first();
    vector& pi = Qpi.second();

    tensor R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
    pi = pi & R;
    Q = Q & R;

    R = rotationTensorY(0.5*deltaT*pi.y()/momentOfInertia_.yy());
    pi = pi & R;
    Q = Q & R;

    R = rotationTensorZ(deltaT*pi.z()/momentOfInertia_.zz());
    pi = pi & R;
    Q = Q & R;

    R = rotationTensorY(0.5*deltaT*pi.y()/momentOfInertia_.yy());
    pi = pi & R;
    Q = Q & R;

    R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
    pi = pi & R;
    Q = Q & R;

    return Qpi;
}

The rotation function below shows that it is actually
0.5*Rx*0.5*Ry*Rz*0.5*Ry*0.5*Rx

[manoj_nav]

Hi Haze

If the orientation tensor is 0.5*Rx*0.5*Ry*Rz*0.5*Ry*0.5*Rx , then the 7th member of the tensor becomes zero for rotation about only y-axis. But in my simulations, I am getting this number non-zero. Can you confirm?



Thanks and Regards,

Manoj

[tmatrix]

Do you know the orientation adjustment from Q0() to Q()?

Originally Posted by haze_1986

For interDyMFoam, if you are using sixDoFRigidBodyMotionI.H:

inline Foam::tensor Foam::sixDoFRigidBodyMotion::rotationTensorX
(
    scalar phi
) const
{
    return tensor
    (
        1, 0, 0,
        0, Foam::cos(phi), -Foam::sin(phi),
        0, Foam::sin(phi), Foam::cos(phi)
    );
}


inline Foam::tensor Foam::sixDoFRigidBodyMotion::rotationTensorY
(
    scalar phi
) const
{
    return tensor
    (
        Foam::cos(phi), 0, Foam::sin(phi),
        0, 1, 0,
        -Foam::sin(phi), 0, Foam::cos(phi)
    );
}


inline Foam::tensor Foam::sixDoFRigidBodyMotion::rotationTensorZ
(
    scalar phi
) const
{
    return tensor
    (
        Foam::cos(phi), -Foam::sin(phi), 0,
        Foam::sin(phi), Foam::cos(phi), 0,
        0, 0, 1
    );
}


inline Foam::Tuple2<Foam::tensor, Foam::vector>
Foam::sixDoFRigidBodyMotion::rotate
(
    const tensor& Q0,
    const vector& pi0,
    const scalar deltaT
) const
{
    Tuple2<tensor, vector> Qpi(Q0, pi0);
    tensor& Q = Qpi.first();
    vector& pi = Qpi.second();

    tensor R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
    pi = pi & R;
    Q = Q & R;

    R = rotationTensorY(0.5*deltaT*pi.y()/momentOfInertia_.yy());
    pi = pi & R;
    Q = Q & R;

    R = rotationTensorZ(deltaT*pi.z()/momentOfInertia_.zz());
    pi = pi & R;
    Q = Q & R;

    R = rotationTensorY(0.5*deltaT*pi.y()/momentOfInertia_.yy());
    pi = pi & R;
    Q = Q & R;

    R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
    pi = pi & R;
    Q = Q & R;

    return Qpi;
}

The rotation function below shows that it is actually
0.5*Rx*0.5*Ry*Rz*0.5*Ry*0.5*Rx

[lil_star]

[QUOTE=manoj_nav;581507]Got this working by including,

awk 'NR == 1 || NR % 3 == 0' centreMass_Extract > centreMass_Extract1


But I still dont understand what are these 9 values of orientation tensor..


where did you add this after sed ?

[zkdkeen]

I think the script should be modified as below:

awk 'NR == 1 || (NR-1) % 3 == 0' centreMass_Extract > centreMass_Extract1

so that the coordinates of centre of mass will match the time correctly.

As to the orientation, I think if you fix your ship in y=0 plane, the orientation should be (cosθ 0 sinθ 0 1 0 −sinθ 0 cosθ).

[lil_star]

In dynamicmeshdict if y-plane = 1 ; then that means the object is allowed to rotate about y-axis. it will have pitch. right?



constraints
{
zAxis
{
sixDoFRigidBodyMotionConstraint line;
direction (0 0 1);
}
yPlane
{
sixDoFRigidBodyMotionConstraint axis;
axis (0 1 0);
}

[mahsankhan]

Originally Posted by lil_star

In dynamicmeshdict if y-plane = 1 ; then that means the object is allowed to rotate about y-axis. it will have pitch. right?



constraints
{
zAxis
{
sixDoFRigidBodyMotionConstraint line;
direction (0 0 1);
}
yPlane
{
sixDoFRigidBodyMotionConstraint axis;
axis (0 1 0);
}

Yes according to this if x-axis is along the wave travel and z is your vertical axis, then the constraints that you have shown are for linear motion in heave and rotation about y in pitch.

But can you please tell me how can we get the values of the pitch in degrees (or even in radians) from the orientation tensor? There is a total of 9 values, I am not understanding it. Any help would be nice.

[mAlletto]

Originally Posted by haze_1986

For interDyMFoam, if you are using sixDoFRigidBodyMotionI.H:

inline Foam::tensor Foam::sixDoFRigidBodyMotion::rotationTensorX
(
    scalar phi
) const
{
    return tensor
    (
        1, 0, 0,
        0, Foam::cos(phi), -Foam::sin(phi),
        0, Foam::sin(phi), Foam::cos(phi)
    );
}


inline Foam::tensor Foam::sixDoFRigidBodyMotion::rotationTensorY
(
    scalar phi
) const
{
    return tensor
    (
        Foam::cos(phi), 0, Foam::sin(phi),
        0, 1, 0,
        -Foam::sin(phi), 0, Foam::cos(phi)
    );
}


inline Foam::tensor Foam::sixDoFRigidBodyMotion::rotationTensorZ
(
    scalar phi
) const
{
    return tensor
    (
        Foam::cos(phi), -Foam::sin(phi), 0,
        Foam::sin(phi), Foam::cos(phi), 0,
        0, 0, 1
    );
}


inline Foam::Tuple2<Foam::tensor, Foam::vector>
Foam::sixDoFRigidBodyMotion::rotate
(
    const tensor& Q0,
    const vector& pi0,
    const scalar deltaT
) const
{
    Tuple2<tensor, vector> Qpi(Q0, pi0);
    tensor& Q = Qpi.first();
    vector& pi = Qpi.second();

    tensor R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
    pi = pi & R;
    Q = Q & R;

    R = rotationTensorY(0.5*deltaT*pi.y()/momentOfInertia_.yy());
    pi = pi & R;
    Q = Q & R;

    R = rotationTensorZ(deltaT*pi.z()/momentOfInertia_.zz());
    pi = pi & R;
    Q = Q & R;

    R = rotationTensorY(0.5*deltaT*pi.y()/momentOfInertia_.yy());
    pi = pi & R;
    Q = Q & R;

    R = rotationTensorX(0.5*deltaT*pi.x()/momentOfInertia_.xx());
    pi = pi & R;
    Q = Q & R;

    return Qpi;
}

The rotation function below shows that it is actually
0.5*Rx*0.5*Ry*Rz*0.5*Ry*0.5*Rx

Does someone know why the rotation Matrix is constructed like this? (0.5*Rx*0.5*Ry*Rz*0.5*Ry*0.5*Rx)


Best


Michael

[umutcaninal]

could you find a solution to that problem?
I am also dealing with the same issue. I have an orientation tensor in the output file as;

(1 0 8.11758850564-22
0 1 0
-8.11758850564-22 0 1)
It's the Ry matrix that gives the pitch angle of the rigid body as it introduced in the link below;
http://planning.cs.uiuc.edu/node102.html

Please let me know if i am mistaken, best regards.

[zyfsoton]

Dear all,

I have the same question about this pitch and heave plot.

I think the sin(pitch) is the third component of the orientation vector but I have no idea how to derive the heave.

Please let me know if you have any idea.

Many thanks,
Tony

[JNSN]

since some years there is a function object available, which writes the motion state to file. No need for parsing log files and doing transformations on your own.


Best,
Jan

[Felipe.td]

Originally Posted by JNSN

since some years there is a function object available, which writes the motion state to file. No need for parsing log files and doing transformations on your own.


Best,
Jan

Can you tell which one It is? And How it Works?
Thanks

[JNSN]

If you are using the rigidBodyDynamics library, then put folowing coding in the functions section of your controlDict:

  rigidBodyState
  {
    angleFormat degrees;
    functionObjectLibs
      (
        "librigidBodyState.so"
      );
    type rigidBodyState;
   }

For the sixDoFRigidBodyMotions lib, adjust the names accordingly.


Best,
Jan