[quarkz]

Hi,

I'm trying to understand how OF does its 6DOF calculation.

In most aircraft 6DOF which uses the equation of motion, they use the body fixed coordinate sys, so the axes changes direction together with the body.

In OF, the forces and moments calculated are fixed in the inertial reference frame, so the axes do not change direction.

In this case, are the 6DOF eqn of motions in OF different from the aircraft 6DOF?

For e.g. Fx = m(u_dot + q*w - r*v)

where u_dot = acc along x dir
q = angular vel in pitch
w = vel in z dir
r = ang. vel in yaw
v = vel in y dir

The additional q*w and r*v terms are due to the body fixed frame.

So in OF, do these additional terms exist? Since OF uses inertial frame.

Hope someone can clarify it.

[Saleh Abuhanieh]

Hi,

As far I remember, in OpenFOAM code (and in many literature), the linear displacements are found using the inertial frame. However, for the angular displacements, it uses the body-fixed frame.

Regards,
Saleh

[quarkz]

Thanks Saleh.

So the additional q*w and r*v terms are not included in the calculation, is that so?

I also looked into the source code:

void Foam::sixDoFRigidBodyMotion::updateAcceleration
(
    const vector& fGlobal,
    const vector& tauGlobal
)
{
    static bool first = true;

    // Save the previous iteration accelerations for relaxation
    vector aPrevIter = a();
    vector tauPrevIter = tau();

    // Calculate new accelerations
    a() = fGlobal/mass_;
    tau() = (Q().T() & tauGlobal);
    applyRestraints();

As shown, a() = fGlobal/mass_

Also, fGlobal represents forces in the global inertial reference frame.

[Saleh Abuhanieh]

Hi,


Yes, the fGlobal includes only the inertial forces.
I've noted down the 6DOF code sequence execution in the past, I am writing it hereunder, hope you can find it useful.



[dynamicMeshDict]

solver sixDoFRigidBodyMotion;

//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//

[sixDoFRigidBodyMotionSolver.C]

Foam::sixDoFRigidBodyMotionSolver::sixDoFRigidBody MotionSolver
(
const polyMesh& mesh,
const IOdictionary& dict
)
:
displacementMotionSolver(mesh, dict, typeName),
motion_ //is a "sixDoFRigidBodyMotion" object;
(
coeffDict(),
IOobject
(
"sixDoFRigidBodyMotionState",
mesh.time().timeName(),
"uniform",
mesh
).typeHeaderOk<IOdictionary>(true)
? IOdictionary
(
IOobject
(
"sixDoFRigidBodyMotionState",
mesh.time().timeName(),
"uniform",
mesh,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE,
false
)
)
: coeffDict(),
mesh.time()
),
patches_(coeffDict().get<wordRes>("patches")),
..



void Foam::sixDoFRigidBodyMotionSolver::solve()
{
..
functionObjects::forces f("forces", db(), forcesDict);

f.calcForcesMoment();

motion_.update
(
firstIter,
ramp*(f.forceEff() + motion_.mass()*g.value()),
ramp
*(
f.momentEff()
+ motion_.mass()*(motion_.momentArm() ^ g.value())
),
t.deltaTValue(),
t.deltaT0Value()
);
..
}


//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//

[sixDoFRigidBodyMotion.C]

void Foam::sixDoFRigidBodyMotion::update
(
bool firstIter,
const vector& fGlobal,
const vector& tauGlobal,
scalar deltaT,
scalar deltaT0
)
{
if (Pstream::master())
{
solver_->solve(firstIter, fGlobal, tauGlobal, deltaT, deltaT0); //integrator, e.g, Newmark

if (report_)
{
status();
}
}

Pstream::scatter(motionState_);
}


//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//

[Newmark.C]

void Foam::sixDoFSolvers::Newmark::solve
(
bool firstIter,
const vector& fGlobal,
const vector& tauGlobal,
scalar deltaT,
scalar deltaT0
)
{
// Update the linear acceleration and torque
updateAcceleration(fGlobal, tauGlobal);

// Correct linear velocity
v() =
tConstraints()
& (v0() + aDamp()*deltaT*(gamma_*a() + (1 - gamma_)*a0()));

// Correct angular momentum
pi() =
rConstraints()
& (pi0() + aDamp()*deltaT*(gamma_*tau() + (1 - gamma_)*tau0()));

// Correct position
centreOfRotation() =
centreOfRotation0()
+ (
tConstraints()
& (
deltaT*v0()
+ aDamp()*sqr(deltaT)*(beta_*a() + (0.5 - beta_)*a0())
)
);

// Correct orientation
vector piDeltaT =
rConstraints()
& (
deltaT*pi0()
+ aDamp()*sqr(deltaT)*(beta_*tau() + (0.5 - beta_)*tau0())
);
Tuple2<tensor, vector> Qpi = rotate(Q0(), piDeltaT, 1);
Q() = Qpi.first(); //update the orientation tensor
}

//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//

[sixDoFSolverI.H]

inline void Foam::sixDoFSolver::updateAcceleration
(
const vector& fGlobal,
const vector& tauGlobal
)
{
body_.updateAcceleration(fGlobal, tauGlobal);
}

//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//

[sixDoFRigidBodyMotion.C]

void Foam::sixDoFRigidBodyMotion::updateAcceleration
(
const vector& fGlobal,
const vector& tauGlobal
)
{
static bool first = true;

// Save the previous iteration accelerations for relaxation
vector aPrevIter = a();
vector tauPrevIter = tau();

// Calculate new accelerations
a() = fGlobal/mass_;
tau() = (Q().T() & tauGlobal);
applyRestraints();

// Relax accelerations on all but first iteration
if (!first)
{
a() = aRelax_*a() + (1 - aRelax_)*aPrevIter;
tau() = aRelax_*tau() + (1 - aRelax_)*tauPrevIter;
}
else
{
first = false;
}
}

//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//

[sixDoFRigidBodyMotionI.H]

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;
}

[quarkz]

Hi Saleh,

Thanks for the info. This is really useful!

The way the 6DOF code sequence works is really confusing. This layout is much easier to understand.