First version of geometricSixDOF class
Compiles but does not produce correct results. Still need to find bugs.
This commit is contained in:
Vuko Vukcevic
Hrvoje Jasak
parent
60821402da
commit
546cfc0f5f
6 changed files with 1107 additions and 44 deletions
|
|
@ -17,6 +17,7 @@ $(sixDOF)/sixDOFODE/sixDOFODEIO.C
|
|||
$(sixDOF)/sixDOFODE/sixDOFODE.C
|
||||
$(sixDOF)/sixDOFODE/newSixDOFODE.C
|
||||
$(sixDOF)/quaternionSixDOF/quaternionSixDOF.C
|
||||
$(sixDOF)/geometricSixDOF/geometricSixDOF.C
|
||||
|
||||
$(sixDOF)/sixDOFBodies/sixDOFBodies.C
|
||||
|
||||
|
|
|
|||
703
src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C
Normal file
703
src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C
Normal file
|
|
@ -0,0 +1,703 @@
|
|||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | foam-extend: Open Source CFD
|
||||
\\ / O peration | Version: 4.0
|
||||
\\ / A nd | Web: http://www.foam-extend.org
|
||||
\\/ M anipulation | For copyright notice see file Copyright
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of foam-extend.
|
||||
|
||||
foam-extend is free software: you can redistribute it and/or modify it
|
||||
under the terms of the GNU General Public License as published by the
|
||||
Free Software Foundation, either version 3 of the License, or (at your
|
||||
option) any later version.
|
||||
|
||||
foam-extend is distributed in the hope that it will be useful, but
|
||||
WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Class
|
||||
geometricSixDOF
|
||||
|
||||
Description
|
||||
6-DOF solver using a geometric method for integration of rotations.
|
||||
|
||||
Author
|
||||
Viktor Pandza, FSB Zagreb. All rights reserved.
|
||||
Vuko Vukcevic, FSB Zagreb. All rights reserved.
|
||||
|
||||
SourceFiles
|
||||
geometricSixDOF.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "geometricSixDOF.H"
|
||||
#include "OutputControlDictionary.H"
|
||||
#include "addToRunTimeSelectionTable.H"
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
defineTypeNameAndDebug(geometricSixDOF, 0);
|
||||
addToRunTimeSelectionTable(sixDOFODE, geometricSixDOF, dictionary);
|
||||
|
||||
}
|
||||
|
||||
|
||||
const Foam::debug::tolerancesSwitch Foam::geometricSixDOF::rotIncTensorTol_
|
||||
(
|
||||
"geometrixSixDOFRotIncTensorTol",
|
||||
1e-10
|
||||
);
|
||||
|
||||
|
||||
const Foam::debug::tolerancesSwitch Foam::geometricSixDOF::rotIncRateTol_
|
||||
(
|
||||
"geometrixSixDOFRotIncRateTol",
|
||||
1e-6
|
||||
);
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
Foam::dimensionedVector Foam::geometricSixDOF::A
|
||||
(
|
||||
const dimensionedVector& xR,
|
||||
const dimensionedVector& uR,
|
||||
const tensor& R
|
||||
) const
|
||||
{
|
||||
// Fix the total force in global coordinate system
|
||||
dimensionedVector fAbs =
|
||||
// Force in global coordinate system
|
||||
force()
|
||||
// Force in local coordinate system
|
||||
+ (dimensionedTensor("R_T", dimless, R.T()) & forceRelative())
|
||||
// Spring force in global coordinate system
|
||||
- (linSpringCoeffs() & xR)
|
||||
// Damping force in global coordinate system
|
||||
- (linDampingCoeffs() & uR);
|
||||
|
||||
// Constrain translation simply by setting the total force to zero
|
||||
constrainTranslation(fAbs.value());
|
||||
|
||||
return fAbs/mass();
|
||||
}
|
||||
|
||||
|
||||
Foam::dimensionedVector Foam::geometricSixDOF::OmegaDot
|
||||
(
|
||||
const tensor& R,
|
||||
const dimensionedVector& omega
|
||||
) const
|
||||
{
|
||||
// External moment (torque) in local coordinate system
|
||||
dimensionedVector mRel =
|
||||
// Moment in global coordinate system
|
||||
(dimensionedTensor("R", dimless, R) & moment())
|
||||
// Moment in local coordinate system
|
||||
+ momentRelative();
|
||||
|
||||
// Note: constraints not implemented at the moment. They shall be
|
||||
// implemented in terms of Lagrange multipliers.
|
||||
|
||||
return
|
||||
inv(momentOfInertia())
|
||||
& (
|
||||
E(omega)
|
||||
+ mRel
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
Foam::dimensionedVector Foam::geometricSixDOF::E
|
||||
(
|
||||
const dimensionedVector& omega
|
||||
) const
|
||||
{
|
||||
return (*(momentOfInertia() & omega) & omega);
|
||||
}
|
||||
|
||||
|
||||
void Foam::geometricSixDOF::constrainTranslation(vector& vec) const
|
||||
{
|
||||
// Constrain the vector with respect to referent or global coordinate system
|
||||
if (referentMotionConstraints_)
|
||||
{
|
||||
vector consVec(referentRotation_.R() & vec);
|
||||
|
||||
if (fixedSurge_)
|
||||
{
|
||||
consVec.x() = 0;
|
||||
}
|
||||
if (fixedSway_)
|
||||
{
|
||||
consVec.y() = 0;
|
||||
}
|
||||
if (fixedHeave_)
|
||||
{
|
||||
consVec.z() = 0;
|
||||
}
|
||||
|
||||
vec = referentRotation_.invR() & consVec;
|
||||
}
|
||||
else
|
||||
{
|
||||
if (fixedSurge_)
|
||||
{
|
||||
vec.x() = 0;
|
||||
}
|
||||
if (fixedSway_)
|
||||
{
|
||||
vec.y() = 0;
|
||||
}
|
||||
if (fixedHeave_)
|
||||
{
|
||||
vec.z() = 0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
Foam::tensor Foam::geometricSixDOF::expMap(const vector& rotInc) const
|
||||
{
|
||||
tensor R;
|
||||
|
||||
// Calculate the magnitude of the rotation increment vector
|
||||
const scalar magRotInc = mag(rotInc);
|
||||
|
||||
if (magRotInc < rotIncTensorTol_)
|
||||
{
|
||||
// No rotational increment
|
||||
R = I;
|
||||
}
|
||||
else
|
||||
{
|
||||
// Calculate rotational increment tensor using the exponential map
|
||||
|
||||
// Skew-symmetric tensor corresponding to the rotation increment
|
||||
const tensor skewRotInc(*rotInc);
|
||||
|
||||
R = I
|
||||
+ skewRotInc*sin(magRotInc)/magRotInc
|
||||
+ (skewRotInc & skewRotInc)*(1.0 - cos(magRotInc))/sqr(magRotInc);
|
||||
}
|
||||
|
||||
return R;
|
||||
}
|
||||
|
||||
|
||||
Foam::vector Foam::geometricSixDOF::dexpMap
|
||||
(
|
||||
const vector& rotInc,
|
||||
const vector& omega
|
||||
) const
|
||||
{
|
||||
vector rotIncDot;
|
||||
|
||||
// Calculate the magnitude of the rotation increment vector
|
||||
const scalar magRotInc = mag(rotInc);
|
||||
|
||||
if (magRotInc < rotIncRateTol_)
|
||||
{
|
||||
// Stabilised calculation of rotation increment to avoid small
|
||||
// denominators
|
||||
rotIncDot = omega;
|
||||
}
|
||||
else
|
||||
{
|
||||
// Calculate rate of the rotational increment vector using the
|
||||
// differential of the exponential map
|
||||
|
||||
// Skew-symmetric tensor corresponding to the rotation increment
|
||||
const tensor skewRotInc(*rotInc);
|
||||
|
||||
rotIncDot =
|
||||
(
|
||||
I
|
||||
+ 0.5*skewRotInc
|
||||
- (skewRotInc & skewRotInc)*
|
||||
(magRotInc/tan(magRotInc/2.0) - 2.0)/(2.0*sqr(magRotInc))
|
||||
)
|
||||
& omega;
|
||||
}
|
||||
|
||||
return rotIncDot;
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::geometricSixDOF::geometricSixDOF(const IOobject& io)
|
||||
:
|
||||
sixDOFODE(io),
|
||||
|
||||
Xrel_(dict().lookup("Xrel")),
|
||||
U_(dict().lookup("U")),
|
||||
Uaverage_(U_),
|
||||
rotation_(tensor(dict().lookup("rotationTensor"))),
|
||||
omega_(dict().lookup("omega")),
|
||||
omegaAverage_(omega_),
|
||||
omegaAverageAbsolute_(omega_),
|
||||
|
||||
nEqns_(),
|
||||
coeffs_(),
|
||||
|
||||
fixedSurge_(dict().lookup("fixedSurge")),
|
||||
fixedSway_(dict().lookup("fixedSway")),
|
||||
fixedHeave_(dict().lookup("fixedHeave")),
|
||||
fixedRoll_(dict().lookup("fixedRoll")),
|
||||
fixedPitch_(dict().lookup("fixedPitch")),
|
||||
fixedYaw_(dict().lookup("fixedYaw")),
|
||||
referentMotionConstraints_
|
||||
(
|
||||
dict().lookupOrDefault<Switch>
|
||||
(
|
||||
"referentMotionConstraints",
|
||||
false
|
||||
)
|
||||
),
|
||||
referentRotation_(vector(1, 0, 0), 0)
|
||||
{
|
||||
// Missing constraints. Count how many rotational constraints we have in
|
||||
// order to count the number of equations.
|
||||
|
||||
nEqns_ = 12;
|
||||
|
||||
// Set size for ODE coefficients depending on number of equations
|
||||
coeffs_.setSize(nEqns_);
|
||||
|
||||
// Set ODE coefficients from position and rotation
|
||||
|
||||
// Linear displacement relative to spring equilibrium
|
||||
const vector& Xval = Xrel_.value();
|
||||
coeffs_[0] = Xval.x();
|
||||
coeffs_[1] = Xval.y();
|
||||
coeffs_[2] = Xval.z();
|
||||
|
||||
// Linear velocity
|
||||
const vector& Uval = U_.value();
|
||||
coeffs_[3] = Uval.x();
|
||||
coeffs_[4] = Uval.y();
|
||||
coeffs_[5] = Uval.z();
|
||||
|
||||
// Rotational velocity in non - inertial coordinate system
|
||||
const vector& omegaVal = omega_.value();
|
||||
coeffs_[6] = omegaVal.x();
|
||||
coeffs_[7] = omegaVal.y();
|
||||
coeffs_[8] = omegaVal.z();
|
||||
|
||||
// Increment of the rotation vector (zero for initial condition)
|
||||
coeffs_[9] = 0;
|
||||
coeffs_[10] = 0;
|
||||
coeffs_[11] = 0;
|
||||
}
|
||||
|
||||
|
||||
Foam::geometricSixDOF::geometricSixDOF
|
||||
(
|
||||
const word& name,
|
||||
const geometricSixDOF& gsd
|
||||
)
|
||||
:
|
||||
sixDOFODE
|
||||
(
|
||||
IOobject
|
||||
(
|
||||
name,
|
||||
gsd.dict().instance(),
|
||||
gsd.dict().local(),
|
||||
gsd.dict().db(),
|
||||
IOobject::NO_READ,
|
||||
IOobject::NO_WRITE
|
||||
)
|
||||
),
|
||||
|
||||
Xrel_(gsd.Xrel_.name(), gsd.Xrel_),
|
||||
U_(gsd.U_.name(), gsd.U_),
|
||||
Uaverage_(gsd.Uaverage_.name(), gsd.Uaverage_),
|
||||
rotation_(gsd.rotation_),
|
||||
omega_(gsd.omega_.name(), gsd.omega_),
|
||||
omegaAverage_(gsd.omegaAverage_.name(), gsd.omegaAverage_),
|
||||
omegaAverageAbsolute_
|
||||
(
|
||||
gsd.omegaAverageAbsolute_.name(),
|
||||
gsd.omegaAverageAbsolute_
|
||||
),
|
||||
|
||||
nEqns_(gsd.nEqns_),
|
||||
coeffs_(gsd.coeffs_),
|
||||
|
||||
fixedSurge_(gsd.fixedSurge_),
|
||||
fixedSway_(gsd.fixedSway_),
|
||||
fixedHeave_(gsd.fixedHeave_),
|
||||
fixedRoll_(gsd.fixedRoll_),
|
||||
fixedPitch_(gsd.fixedPitch_),
|
||||
fixedYaw_(gsd.fixedYaw_),
|
||||
referentMotionConstraints_(gsd.referentMotionConstraints_),
|
||||
referentRotation_(gsd.referentRotation_)
|
||||
{}
|
||||
|
||||
|
||||
Foam::autoPtr<Foam::sixDOFODE> Foam::geometricSixDOF::clone
|
||||
(
|
||||
const word& name
|
||||
) const
|
||||
{
|
||||
// Create and return an autoPtr to the new geometricSixDOF object with a
|
||||
// new name
|
||||
return autoPtr<sixDOFODE>
|
||||
(
|
||||
new geometricSixDOF
|
||||
(
|
||||
name,
|
||||
*this
|
||||
)
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::geometricSixDOF::~geometricSixDOF()
|
||||
{}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
||||
|
||||
const Foam::dimensionedVector& Foam::geometricSixDOF::Xrel() const
|
||||
{
|
||||
return Xrel_;
|
||||
}
|
||||
|
||||
|
||||
const Foam::dimensionedVector& Foam::geometricSixDOF::omega() const
|
||||
{
|
||||
return omega_;
|
||||
}
|
||||
|
||||
|
||||
Foam::dimensionedVector Foam::geometricSixDOF::X() const
|
||||
{
|
||||
return Xequilibrium() + Xrel_;
|
||||
}
|
||||
|
||||
|
||||
const Foam::dimensionedVector& Foam::geometricSixDOF::U() const
|
||||
{
|
||||
return U_;
|
||||
}
|
||||
|
||||
const Foam::dimensionedVector& Foam::geometricSixDOF::Uaverage() const
|
||||
{
|
||||
return Uaverage_;
|
||||
}
|
||||
|
||||
|
||||
const Foam::finiteRotation& Foam::geometricSixDOF::rotation() const
|
||||
{
|
||||
return rotation_;
|
||||
}
|
||||
|
||||
|
||||
Foam::vector Foam::geometricSixDOF::rotVector() const
|
||||
{
|
||||
return rotation_.rotVector();
|
||||
}
|
||||
|
||||
|
||||
Foam::dimensionedScalar Foam::geometricSixDOF::rotAngle() const
|
||||
{
|
||||
return dimensionedScalar("rotAngle", dimless, rotation_.rotAngle());
|
||||
}
|
||||
|
||||
|
||||
void Foam::geometricSixDOF::setXrel(const vector& x)
|
||||
{
|
||||
Xrel_.value() = x;
|
||||
|
||||
// Set X in coefficients
|
||||
coeffs_[0] = x.x();
|
||||
coeffs_[1] = x.y();
|
||||
coeffs_[2] = x.z();
|
||||
}
|
||||
|
||||
|
||||
void Foam::geometricSixDOF::setU(const vector& U)
|
||||
{
|
||||
U_.value() = U;
|
||||
Uaverage_ = U_;
|
||||
|
||||
// Set U in coefficients
|
||||
coeffs_[3] = U.x();
|
||||
coeffs_[4] = U.y();
|
||||
coeffs_[5] = U.z();
|
||||
}
|
||||
|
||||
|
||||
void Foam::geometricSixDOF::setRotation(const HamiltonRodriguezRot& rot)
|
||||
{
|
||||
// Set increment rotation vector to zero
|
||||
coeffs_[9] = 0;
|
||||
coeffs_[10] = 0;
|
||||
coeffs_[11] = 0;
|
||||
}
|
||||
|
||||
|
||||
void Foam::geometricSixDOF::setOmega(const vector& omega)
|
||||
{
|
||||
omega_.value() = omega;
|
||||
omegaAverage_ = omega_;
|
||||
omegaAverageAbsolute_ = omega_;
|
||||
|
||||
// Set omega in coefficients
|
||||
coeffs_[6] = omega.x();
|
||||
coeffs_[7] = omega.y();
|
||||
coeffs_[8] = omega.z();
|
||||
}
|
||||
|
||||
|
||||
void Foam::geometricSixDOF::setReferentRotation
|
||||
(
|
||||
const HamiltonRodriguezRot& rot
|
||||
)
|
||||
{
|
||||
referentRotation_ = rot;
|
||||
|
||||
referentMotionConstraints_ = true;
|
||||
}
|
||||
|
||||
|
||||
void Foam::geometricSixDOF::setState(const sixDOFODE& sd)
|
||||
{
|
||||
// First set the state in base class subobject
|
||||
sixDOFODE::setState(sd);
|
||||
|
||||
// Cast sixDOFODE& to geometricSixDOF&
|
||||
const geometricSixDOF& gsd = refCast<const geometricSixDOF>(sd);
|
||||
|
||||
// Set state variables for this class
|
||||
Xrel_ = gsd.Xrel_;
|
||||
U_ = gsd.U_;
|
||||
Uaverage_ = gsd.Uaverage_;
|
||||
rotation_ = gsd.rotation_;
|
||||
omega_ = gsd.omega_;
|
||||
omegaAverage_ = gsd.omegaAverage_;
|
||||
omegaAverageAbsolute_ = gsd.omegaAverageAbsolute_;
|
||||
|
||||
// Copy ODE coefficients: this carries actual ODE state
|
||||
// HJ, 23/Mar/2015
|
||||
coeffs_ = gsd.coeffs_;
|
||||
|
||||
fixedSurge_ = gsd.fixedSurge_;
|
||||
fixedSway_ = gsd.fixedSway_;
|
||||
fixedHeave_ = gsd.fixedHeave_;
|
||||
fixedRoll_ = gsd.fixedRoll_;
|
||||
fixedPitch_ = gsd.fixedPitch_;
|
||||
fixedYaw_ = gsd.fixedYaw_;
|
||||
referentMotionConstraints_ = gsd.referentMotionConstraints_;
|
||||
referentRotation_ = gsd.referentRotation_;
|
||||
}
|
||||
|
||||
|
||||
Foam::vector Foam::geometricSixDOF::rotVectorAverage() const
|
||||
{
|
||||
return rotation_.rotVectorAverage();
|
||||
}
|
||||
|
||||
|
||||
const Foam::dimensionedVector& Foam::geometricSixDOF::omegaAverage() const
|
||||
{
|
||||
return omegaAverage_;
|
||||
}
|
||||
|
||||
|
||||
const Foam::dimensionedVector&
|
||||
Foam::geometricSixDOF::omegaAverageAbsolute() const
|
||||
{
|
||||
return omegaAverageAbsolute_;
|
||||
}
|
||||
|
||||
|
||||
Foam::tensor Foam::geometricSixDOF::toRelative() const
|
||||
{
|
||||
// Note: using rotation tensor directly since we are integrating rotational
|
||||
// increment vector
|
||||
return rotation_.rotTensor();
|
||||
}
|
||||
|
||||
|
||||
Foam::tensor Foam::geometricSixDOF::toAbsolute() const
|
||||
{
|
||||
// Note: using rotation tensor directly since we are integrating rotational
|
||||
// increment vector
|
||||
return rotation_.rotTensor().T();
|
||||
}
|
||||
|
||||
|
||||
const Foam::tensor& Foam::geometricSixDOF::rotIncrementTensor() const
|
||||
{
|
||||
return rotation_.rotIncrementTensor();
|
||||
}
|
||||
|
||||
|
||||
void Foam::geometricSixDOF::derivatives
|
||||
(
|
||||
const scalar x,
|
||||
const scalarField& y,
|
||||
scalarField& dydx
|
||||
) const
|
||||
{
|
||||
// Translation
|
||||
|
||||
// Set the derivatives for displacement
|
||||
dydx[0] = y[3];
|
||||
dydx[1] = y[4];
|
||||
dydx[2] = y[5];
|
||||
|
||||
dimensionedVector curX("curX", dimLength, vector(y[0], y[1], y[2]));
|
||||
dimensionedVector curU("curU", dimVelocity, vector(y[3], y[4], y[5]));
|
||||
|
||||
// Get rotational increment vector (u)
|
||||
const vector rotIncrementVector(y[9], y[10], y[11]);
|
||||
|
||||
// Calculate rotation increment tensor obtained with exponential map using
|
||||
// the increment vector
|
||||
const tensor rotIncrement = expMap(rotIncrementVector);
|
||||
|
||||
// Update rotation tensor
|
||||
rotation_.updateRotation(rotIncrement);
|
||||
|
||||
// Calculate translational acceleration using current rotation
|
||||
const vector accel = A(curX, curU, rotation_.rotTensor()).value();
|
||||
|
||||
// Set the derivatives for velocity
|
||||
dydx[3] = accel.x();
|
||||
dydx[4] = accel.y();
|
||||
dydx[5] = accel.z();
|
||||
|
||||
// Rotation
|
||||
|
||||
dimensionedVector curOmega
|
||||
(
|
||||
"curOmega",
|
||||
dimless/dimTime,
|
||||
vector(y[6], y[7], y[8])
|
||||
);
|
||||
|
||||
// Calculate rotational acceleration using current rotation
|
||||
const vector omegaDot = OmegaDot(rotation_.rotTensor(), curOmega).value();
|
||||
|
||||
dydx[6] = omegaDot.x();
|
||||
dydx[7] = omegaDot.y();
|
||||
dydx[8] = omegaDot.z();
|
||||
|
||||
// Calculate derivative of rotIncrementVector using current rotation
|
||||
// information
|
||||
const vector rotIncrementVectorDot =
|
||||
dexpMap(rotIncrementVector, curOmega.value());
|
||||
|
||||
// Set the derivatives for rotation
|
||||
dydx[9] = rotIncrementVectorDot.x();
|
||||
dydx[10] = rotIncrementVectorDot.y();
|
||||
dydx[11] = rotIncrementVectorDot.z();
|
||||
}
|
||||
|
||||
|
||||
void Foam::geometricSixDOF::update(const scalar delta)
|
||||
{
|
||||
// Translation
|
||||
|
||||
// Update displacement
|
||||
vector Xold = Xrel_.value();
|
||||
|
||||
vector& Xval = Xrel_.value();
|
||||
|
||||
Xval.x() = coeffs_[0];
|
||||
Xval.y() = coeffs_[1];
|
||||
Xval.z() = coeffs_[2];
|
||||
|
||||
// Update velocity
|
||||
Uaverage_.value() = (Xval - Xold)/delta;
|
||||
|
||||
vector& Uval = U_.value();
|
||||
|
||||
Uval.x() = coeffs_[3];
|
||||
Uval.y() = coeffs_[4];
|
||||
Uval.z() = coeffs_[5];
|
||||
|
||||
// Constrain velocity and re-set coefficients
|
||||
constrainTranslation(Uval);
|
||||
coeffs_[3] = Uval.x();
|
||||
coeffs_[4] = Uval.y();
|
||||
coeffs_[5] = Uval.z();
|
||||
|
||||
// Rotation
|
||||
|
||||
// Update omega
|
||||
vector& omegaVal = omega_.value();
|
||||
|
||||
omegaVal.x() = coeffs_[6];
|
||||
omegaVal.y() = coeffs_[7];
|
||||
omegaVal.z() = coeffs_[8];
|
||||
|
||||
// Update rotation with final increment vector
|
||||
rotation_.updateRotation
|
||||
(
|
||||
expMap(vector(coeffs_[9], coeffs_[10], coeffs_[11]))
|
||||
);
|
||||
|
||||
// Reset increment vector in coefficients for the next step
|
||||
coeffs_[9] = 0;
|
||||
coeffs_[10] = 0;
|
||||
coeffs_[11] = 0;
|
||||
|
||||
omegaAverage_.value() = rotation_.omegaAverage(delta);
|
||||
|
||||
// Calculate omegaAverageAbsolute
|
||||
omegaAverageAbsolute_.value() = rotation_.omegaAverageAbsolute(delta);
|
||||
}
|
||||
|
||||
|
||||
bool Foam::geometricSixDOF::writeData(Ostream& os) const
|
||||
{
|
||||
// First write the part related to base class subobject
|
||||
sixDOFODE::writeData(os);
|
||||
|
||||
// Write type name
|
||||
os.writeKeyword("type") << tab << type() << token::END_STATEMENT << endl;
|
||||
|
||||
// Write data
|
||||
os.writeKeyword("Xrel") << tab << this->Xrel()
|
||||
<< token::END_STATEMENT << nl;
|
||||
os.writeKeyword("U") << tab << this->U() << token::END_STATEMENT << nl;
|
||||
os.writeKeyword("rotationTensor") << tab << this->toRelative()
|
||||
<< token::END_STATEMENT << nl;
|
||||
os.writeKeyword("omega") << tab << this->omega()
|
||||
<< token::END_STATEMENT << nl << nl;
|
||||
|
||||
os.writeKeyword("fixedSurge") << tab << this->fixedSurge_ <<
|
||||
token::END_STATEMENT << endl;
|
||||
os.writeKeyword("fixedSway") << tab << this->fixedSway_ <<
|
||||
token::END_STATEMENT << endl;
|
||||
os.writeKeyword("fixedHeave") << tab << this->fixedHeave_ <<
|
||||
token::END_STATEMENT << endl;
|
||||
os.writeKeyword("fixedRoll") << tab << this->fixedRoll_ <<
|
||||
token::END_STATEMENT << endl;
|
||||
os.writeKeyword("fixedPitch") << tab << this->fixedPitch_ <<
|
||||
token::END_STATEMENT << endl;
|
||||
os.writeKeyword("fixedYaw") << tab << this->fixedYaw_ <<
|
||||
token::END_STATEMENT << endl;
|
||||
|
||||
return os.good();
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
378
src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H
Normal file
378
src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H
Normal file
|
|
@ -0,0 +1,378 @@
|
|||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / F ield | foam-extend: Open Source CFD
|
||||
\\ / O peration | Version: 4.0
|
||||
\\ / A nd | Web: http://www.foam-extend.org
|
||||
\\/ M anipulation | For copyright notice see file Copyright
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of foam-extend.
|
||||
|
||||
foam-extend is free software: you can redistribute it and/or modify it
|
||||
under the terms of the GNU General Public License as published by the
|
||||
Free Software Foundation, either version 3 of the License, or (at your
|
||||
option) any later version.
|
||||
|
||||
foam-extend is distributed in the hope that it will be useful, but
|
||||
WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
||||
General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with foam-extend. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
Class
|
||||
geometricSixDOF
|
||||
|
||||
Description
|
||||
6-DOF solver using a geometric method for integration of rotations.
|
||||
|
||||
Reference (bibtex):
|
||||
|
||||
@article {mullerTerze2016,
|
||||
Author = {M\"{u}ller, A. and Terze, Z.},
|
||||
title = {Geometric methods and formulations in computational
|
||||
multibody systems},
|
||||
Journal = {Acta Mechanica},
|
||||
Year = {2016},
|
||||
Volume = {227},
|
||||
Number = {12},
|
||||
Pages = {3327--3350}
|
||||
}
|
||||
|
||||
Author
|
||||
Viktor Pandza, FSB Zagreb. All rights reserved.
|
||||
Vuko Vukcevic, FSB Zagreb. All rights reserved.
|
||||
|
||||
SourceFiles
|
||||
geometricSixDOF.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef geometricSixDOF_H
|
||||
#define geometricSixDOF_H
|
||||
|
||||
#include "sixDOFODE.H"
|
||||
#include "finiteRotation.H"
|
||||
#include "tolerancesSwitch.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class geometricSixDOF Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class geometricSixDOF
|
||||
:
|
||||
public sixDOFODE
|
||||
{
|
||||
// Private data
|
||||
|
||||
// Initial body state variables
|
||||
|
||||
//- Displacement relative to spring equilibrium
|
||||
dimensionedVector Xrel_;
|
||||
|
||||
//- Velocity of mass centroid
|
||||
dimensionedVector U_;
|
||||
|
||||
//- Average velocity of mass centroid at previous time-step
|
||||
dimensionedVector Uaverage_;
|
||||
|
||||
//- Finite rotation. Note: changed during solution process
|
||||
mutable finiteRotation rotation_;
|
||||
|
||||
//- Rotational velocity about mass centroid
|
||||
dimensionedVector omega_;
|
||||
|
||||
|
||||
// Average variables that need to be stored
|
||||
|
||||
//- Average rotational velocity in relative coordinate system
|
||||
dimensionedVector omegaAverage_;
|
||||
|
||||
//- Average rotational velocity in absolute coordinate system
|
||||
dimensionedVector omegaAverageAbsolute_;
|
||||
|
||||
|
||||
// ODE controls
|
||||
|
||||
//- Number of equations (depending on rotational constraints)
|
||||
scalar nEqns_;
|
||||
|
||||
//- ODE coefficients
|
||||
scalarField coeffs_;
|
||||
|
||||
|
||||
//- Motion constraints (given as fixed motion components)
|
||||
|
||||
//- Fixed surge (x-translation)
|
||||
Switch fixedSurge_;
|
||||
|
||||
//- Fixed sway (y-translation)
|
||||
Switch fixedSway_;
|
||||
|
||||
//- Fixed heave (z-translation)
|
||||
Switch fixedHeave_;
|
||||
|
||||
//- Fixed roll (rotation around x)
|
||||
Switch fixedRoll_;
|
||||
|
||||
//- Fixed pitch (rotation around y)
|
||||
Switch fixedPitch_;
|
||||
|
||||
//- Fixed yaw (rotation around z)
|
||||
Switch fixedYaw_;
|
||||
|
||||
//- Constraints in referent coordinate system
|
||||
Switch referentMotionConstraints_;
|
||||
|
||||
//- Rotation of referent coordinate system. Currently defined using
|
||||
// quaternions for comatibility. Need to reformulate.
|
||||
// VV, 28/Feb/2017.
|
||||
HamiltonRodriguezRot referentRotation_;
|
||||
|
||||
|
||||
// Private Member Functions
|
||||
|
||||
//- Disallow default bitwise copy construct
|
||||
geometricSixDOF(const geometricSixDOF&);
|
||||
|
||||
//- Disallow default bitwise assignment
|
||||
void operator=(const geometricSixDOF&);
|
||||
|
||||
|
||||
// Variables in relative coordinate system (solved for)
|
||||
|
||||
//- Return acceleration in relative coordinate system
|
||||
// given current values of relative displacement and velocity
|
||||
dimensionedVector A
|
||||
(
|
||||
const dimensionedVector& xR,
|
||||
const dimensionedVector& uR,
|
||||
const tensor& R
|
||||
) const;
|
||||
|
||||
|
||||
//- Return rotational acceleration in relative coordinate system
|
||||
// given current values for relative rotational velocity
|
||||
dimensionedVector OmegaDot
|
||||
(
|
||||
const tensor& R,
|
||||
const dimensionedVector& omega
|
||||
) const;
|
||||
|
||||
//- Return the Euler part of moment equation
|
||||
dimensionedVector E
|
||||
(
|
||||
const dimensionedVector& omega
|
||||
) const;
|
||||
|
||||
//- Constrain translation vector in referent or global coordinate
|
||||
// system
|
||||
void constrainTranslation(vector& vec) const;
|
||||
|
||||
//- Exponential map used to calculate increment of the rotation
|
||||
// tensor
|
||||
tensor expMap(const vector& rotInc) const;
|
||||
|
||||
//- Differential of the expontential map used to calculate the time
|
||||
// derivative of rotation increment vector
|
||||
vector dexpMap(const vector& rotInc, const vector& omega) const;
|
||||
|
||||
|
||||
public:
|
||||
|
||||
// Run-time type information
|
||||
TypeName("geometricSixDOF");
|
||||
|
||||
|
||||
// Static data members
|
||||
|
||||
//- Rotational increment tensor tolerance. Used in expMap member
|
||||
// function in case the rotation is negligibly small
|
||||
static const debug::tolerancesSwitch rotIncTensorTol_;
|
||||
|
||||
//- Rotational increment rate of change tolerance. Used in dexpMap
|
||||
// member function in case the rotation rate is negligibly small
|
||||
static const debug::tolerancesSwitch rotIncRateTol_;
|
||||
|
||||
|
||||
// Constructors
|
||||
|
||||
//- Construct from dictionary
|
||||
geometricSixDOF(const IOobject& io);
|
||||
|
||||
//- Construct geometricSixDOF object, changing name
|
||||
geometricSixDOF
|
||||
(
|
||||
const word& name,
|
||||
const geometricSixDOF& gsd
|
||||
);
|
||||
|
||||
//- Return a clone, changing the name
|
||||
virtual autoPtr<sixDOFODE> clone(const word& name) const;
|
||||
|
||||
|
||||
// Destructor
|
||||
|
||||
virtual ~geometricSixDOF();
|
||||
|
||||
|
||||
// Member Functions
|
||||
|
||||
// Virtual interface for 6DOF motion state
|
||||
|
||||
// Variables in relative coordinate system
|
||||
|
||||
//- Return displacement in translated coordinate system
|
||||
// relative to spring equilibrium
|
||||
virtual const dimensionedVector& Xrel() const;
|
||||
|
||||
//- Return rotational velocity in relative coordinate system
|
||||
virtual const dimensionedVector& omega() const;
|
||||
|
||||
|
||||
// Displacement and rotation in the absolute coordinate system
|
||||
|
||||
//- Return position of origin in absolute coordinate system
|
||||
virtual dimensionedVector X() const;
|
||||
|
||||
//- Return velocity of origin
|
||||
virtual const dimensionedVector& U() const;
|
||||
|
||||
//- Return average velocity of origin for the previous time-step
|
||||
virtual const dimensionedVector& Uaverage() const;
|
||||
|
||||
//- Return finite rotation
|
||||
virtual const finiteRotation& rotation() const;
|
||||
|
||||
//- Return rotational vector of body
|
||||
virtual vector rotVector() const;
|
||||
|
||||
//- Return rotation angle of body
|
||||
virtual dimensionedScalar rotAngle() const;
|
||||
|
||||
|
||||
// Non-access control for setting state variables
|
||||
|
||||
//- Set position of origin
|
||||
virtual void setXrel(const vector& x);
|
||||
|
||||
//- Set velocity of origin
|
||||
virtual void setU(const vector& u);
|
||||
|
||||
//- Set rotational angle in relative coordinate system
|
||||
virtual void setRotation(const HamiltonRodriguezRot& rot);
|
||||
|
||||
//- Set rotational velocity in relative coordinate system
|
||||
virtual void setOmega(const vector& omega);
|
||||
|
||||
//- Set referent coordinate system to apply constraints
|
||||
virtual void setReferentRotation
|
||||
(
|
||||
const HamiltonRodriguezRot& rot
|
||||
);
|
||||
|
||||
//- Set ODE parameters from another ODE
|
||||
virtual void setState(const sixDOFODE&);
|
||||
|
||||
|
||||
// Average motion per time-step
|
||||
|
||||
//- Return average rotational vector of body
|
||||
virtual vector rotVectorAverage() const;
|
||||
|
||||
//- Return average rotational velocity in relative coordinate
|
||||
// system
|
||||
virtual const dimensionedVector& omegaAverage() const;
|
||||
|
||||
//- Return average rotational velocity in absolute coordinate
|
||||
// system
|
||||
virtual const dimensionedVector& omegaAverageAbsolute() const;
|
||||
|
||||
|
||||
// Rotations
|
||||
|
||||
//- Return rotation tensor to relative coordinate system
|
||||
virtual tensor toRelative() const;
|
||||
|
||||
//- Return rotation tensor to absolute coordinate system
|
||||
virtual tensor toAbsolute() const;
|
||||
|
||||
//- Return transformation tensor between new and previous
|
||||
// rotation
|
||||
virtual const tensor& rotIncrementTensor() const;
|
||||
|
||||
|
||||
// ODE parameters
|
||||
|
||||
//- Return number of equations
|
||||
virtual label nEqns() const
|
||||
{
|
||||
return nEqns_;
|
||||
}
|
||||
|
||||
//- Return access to coefficients
|
||||
virtual scalarField& coeffs()
|
||||
{
|
||||
return coeffs_;
|
||||
}
|
||||
|
||||
//- Return reference to coefficients
|
||||
virtual const scalarField& coeffs() const
|
||||
{
|
||||
return coeffs_;
|
||||
}
|
||||
|
||||
//- Evaluate derivatives
|
||||
virtual void derivatives
|
||||
(
|
||||
const scalar x,
|
||||
const scalarField& y,
|
||||
scalarField& dydx
|
||||
) const;
|
||||
|
||||
//- Evaluate Jacobian
|
||||
virtual void jacobian
|
||||
(
|
||||
const scalar x,
|
||||
const scalarField& y,
|
||||
scalarField& dfdx,
|
||||
scalarSquareMatrix& dfdy
|
||||
) const
|
||||
{
|
||||
notImplemented
|
||||
(
|
||||
"geometricSixDOF::jacobian\n"
|
||||
"(\n"
|
||||
" const scalar x,\n"
|
||||
" const scalarField& y,\n"
|
||||
" scalarField& dfdx,\n"
|
||||
" scalarSquareMatrix& dfdy,\n"
|
||||
") const"
|
||||
);
|
||||
}
|
||||
|
||||
//- Update ODE after the solution, advancing by delta
|
||||
virtual void update(const scalar delta);
|
||||
|
||||
// Write controls
|
||||
|
||||
//- WriteData member function required by regIOobject
|
||||
virtual bool writeData(Ostream&) const;
|
||||
};
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
} // End namespace Foam
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
#endif
|
||||
|
||||
// ************************************************************************* //
|
||||
|
|
@ -189,38 +189,6 @@ void Foam::quaternionSixDOF::constrainTranslation(vector& vec) const
|
|||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
|
||||
|
||||
void Foam::quaternionSixDOF::setCoeffs()
|
||||
{
|
||||
// Set ODE coefficients from position and rotation
|
||||
|
||||
// Linear displacement relative to spring equilibrium
|
||||
const vector& Xval = Xrel_.value();
|
||||
coeffs_[0] = Xval.x();
|
||||
coeffs_[1] = Xval.y();
|
||||
coeffs_[2] = Xval.z();
|
||||
|
||||
// Linear velocity
|
||||
const vector& Uval = U_.value();
|
||||
coeffs_[3] = Uval.x();
|
||||
coeffs_[4] = Uval.y();
|
||||
coeffs_[5] = Uval.z();
|
||||
|
||||
// Rotational velocity in non - inertial coordinate system
|
||||
const vector& omegaVal = omega_.value();
|
||||
coeffs_[6] = omegaVal.x();
|
||||
coeffs_[7] = omegaVal.y();
|
||||
coeffs_[8] = omegaVal.z();
|
||||
|
||||
// Quaternions
|
||||
coeffs_[9] = rotation_.eInitial().e0();
|
||||
coeffs_[10] = rotation_.eInitial().e1();
|
||||
coeffs_[11] = rotation_.eInitial().e2();
|
||||
coeffs_[12] = rotation_.eInitial().e3();
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
||||
|
||||
Foam::quaternionSixDOF::quaternionSixDOF(const IOobject& io)
|
||||
|
|
@ -257,7 +225,31 @@ Foam::quaternionSixDOF::quaternionSixDOF(const IOobject& io)
|
|||
),
|
||||
referentRotation_(vector(1, 0, 0), 0)
|
||||
{
|
||||
setCoeffs();
|
||||
// Set ODE coefficients from position and rotation
|
||||
|
||||
// Linear displacement relative to spring equilibrium
|
||||
const vector& Xval = Xrel_.value();
|
||||
coeffs_[0] = Xval.x();
|
||||
coeffs_[1] = Xval.y();
|
||||
coeffs_[2] = Xval.z();
|
||||
|
||||
// Linear velocity
|
||||
const vector& Uval = U_.value();
|
||||
coeffs_[3] = Uval.x();
|
||||
coeffs_[4] = Uval.y();
|
||||
coeffs_[5] = Uval.z();
|
||||
|
||||
// Rotational velocity in non - inertial coordinate system
|
||||
const vector& omegaVal = omega_.value();
|
||||
coeffs_[6] = omegaVal.x();
|
||||
coeffs_[7] = omegaVal.y();
|
||||
coeffs_[8] = omegaVal.z();
|
||||
|
||||
// Quaternions
|
||||
coeffs_[9] = rotation_.eInitial().e0();
|
||||
coeffs_[10] = rotation_.eInitial().e1();
|
||||
coeffs_[11] = rotation_.eInitial().e2();
|
||||
coeffs_[12] = rotation_.eInitial().e3();
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -160,14 +160,6 @@ class quaternionSixDOF
|
|||
void constrainTranslation(vector& vec) const;
|
||||
|
||||
|
||||
protected:
|
||||
|
||||
// Protected Member Functions
|
||||
|
||||
//- Set ODE coefficients from position and rotation
|
||||
virtual void setCoeffs();
|
||||
|
||||
|
||||
public:
|
||||
|
||||
// Run-time type information
|
||||
|
|
|
|||
|
|
@ -146,9 +146,6 @@ protected:
|
|||
|
||||
// Protected Member Functions
|
||||
|
||||
//- Set ODE coefficients from position and rotation
|
||||
virtual void setCoeffs() = 0;
|
||||
|
||||
//- Update Aitkens relaxation parameters
|
||||
void aitkensRelaxation(const scalar min, const scalar max);
|
||||
|
||||
|
|
|
|||
Reference in a new issue