Restructured sixDOFODE interface

Removed implementation specific data and left only essential interface for
coupling with CFD solver.
This commit is contained in:
Vuko Vukcevic 2017-03-01 10:58:59 +01:00
parent 89ef2f50fc
commit d4ed7d23d1
8 changed files with 97 additions and 377 deletions

View file

@ -51,6 +51,9 @@ int main(int argc, char *argv[])
sixDOFBodies structure(runTime);
OFstream of(runTime.path()/"motion.dat");
// Write header for output file
of << "# Time, CoG_x, CoG_y, CoG_z, omega_x, omega_y, omega_z" << endl;
Info<< "\nStarting time loop\n" << endl;
for (runTime++; !runTime.end(); runTime++)
@ -73,11 +76,14 @@ int main(int argc, char *argv[])
<< structure()[bodyI].omegaAverage().value() << nl
<< "Current omega in time instant = "
<< structure()[bodyI].omega().value() << nl
<< "Average omegaAbs in time step = "
<< structure()[bodyI].omegaAverageAbsolute().value() << nl
<< endl;
of << structure()[bodyI].X().value().x() << tab;
of << structure()[bodyI].X().value().x() << tab
<< structure()[bodyI].X().value().y() << tab
<< structure()[bodyI].X().value().z() << tab
<< structure()[bodyI].omega().value().x() << tab
<< structure()[bodyI].omega().value().y() << tab
<< structure()[bodyI].omega().value().z() << tab;
}
of << endl;

View file

@ -40,9 +40,7 @@ Foam::vector Foam::finiteRotation::rotVector(const tensor& rotT)
{
vector ur = - *( inv(I + rotT) & (I - rotT) );
// Scaling to a unit vector. HJ, problems with round-off
// HJ, 4/Aug/2008
// Scaling to a unit vector. Problems with round-off. HJ, 4/Aug/2008
if (mag(ur) > SMALL)
{
return ur/(mag(ur) + SMALL);
@ -50,8 +48,7 @@ Foam::vector Foam::finiteRotation::rotVector(const tensor& rotT)
else
{
// Rotation vector is undertermined at zero rotation
// Returning arbitrary unit vector
// HJ, 4/Mar/2015
// Returning arbitrary unit vector. HJ, 4/Mar/2015
return vector(0, 0, 1);
}
}

View file

@ -170,7 +170,7 @@ Foam::tensor Foam::geometricSixDOF::expMap(const vector& rotInc) const
{
tensor R;
// Calculate the magnitude of the rotation increment vector
// Calculate the magnitude of the rotational increment vector
const scalar magRotInc = mag(rotInc);
if (magRotInc < rotIncTensorTol_)
@ -202,7 +202,7 @@ Foam::vector Foam::geometricSixDOF::dexpMap
{
vector rotIncDot;
// Calculate the magnitude of the rotation increment vector
// Calculate the magnitude of the rotational increment vector
const scalar magRotInc = mag(rotInc);
if (magRotInc < rotIncRateTol_)
@ -241,11 +241,14 @@ Foam::geometricSixDOF::geometricSixDOF(const IOobject& io)
Xrel_(dict().lookup("Xrel")),
U_(dict().lookup("U")),
Uaverage_(U_),
Uaverage_("Uaverage", U_),
rotation_(tensor(dict().lookup("rotationTensor"))),
rotIncrement_
(
dict().lookupOrDefault<tensor>("rotationIncrementTensor", tensor::zero)
),
omega_(dict().lookup("omega")),
omegaAverage_(omega_),
omegaAverageAbsolute_(omega_),
omegaAverage_("omegaAverage", omega_),
nEqns_(),
coeffs_(),
@ -326,11 +329,6 @@ Foam::geometricSixDOF::geometricSixDOF
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_),
@ -401,80 +399,6 @@ const Foam::dimensionedVector& Foam::geometricSixDOF::Uaverage() const
}
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
@ -490,7 +414,6 @@ void Foam::geometricSixDOF::setState(const sixDOFODE& sd)
rotation_ = gsd.rotation_;
omega_ = gsd.omega_;
omegaAverage_ = gsd.omegaAverage_;
omegaAverageAbsolute_ = gsd.omegaAverageAbsolute_;
// Copy ODE coefficients: this carries actual ODE state
// HJ, 23/Mar/2015
@ -507,44 +430,27 @@ void Foam::geometricSixDOF::setState(const sixDOFODE& sd)
}
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();
return rotation_;
}
Foam::tensor Foam::geometricSixDOF::toAbsolute() const
{
// Note: using rotation tensor directly since we are integrating rotational
// increment vector
return rotation_.rotTensor().T();
return rotation_.T();
}
const Foam::tensor& Foam::geometricSixDOF::rotIncrementTensor() const
{
return rotation_.rotIncrementTensor();
return rotIncrement_;
}
@ -569,7 +475,7 @@ void Foam::geometricSixDOF::derivatives
const vector rotIncrementVector(y[9], y[10], y[11]);
// Calculate current rotation tensor obtained with exponential map
const tensor curRot = (expMap(rotIncrementVector) & rotation_.rotTensor());
const tensor curRot = (expMap(rotIncrementVector) & rotation_);
// Calculate translational acceleration using current rotation
const vector accel = A(curX, curU, curRot).value();
@ -612,7 +518,7 @@ void Foam::geometricSixDOF::update(const scalar delta)
// Translation
// Update displacement
vector Xold = Xrel_.value();
const vector Xold = Xrel_.value();
vector& Xval = Xrel_.value();
@ -638,27 +544,28 @@ void Foam::geometricSixDOF::update(const scalar delta)
// Rotation
// Update omega
const vector omegaOld = omega_.value();
vector& omegaVal = omega_.value();
omegaVal.x() = coeffs_[6];
omegaVal.y() = coeffs_[7];
omegaVal.z() = coeffs_[8];
// Update rotation with final increment vector
rotation_.updateRotationWithIncrement
(
expMap(vector(coeffs_[9], coeffs_[10], coeffs_[11]))
);
// Update rotational increment tensor
rotIncrement_ = expMap(vector(coeffs_[9], coeffs_[10], coeffs_[11]));
// Update rotational tensor
rotation_ = (rotIncrement_ & rotation_);
// 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);
// Consider calculating average omega using rotational tensor and rotational
// increment tensors
omegaAverage_.value() = 0.5*(omegaVal + omegaOld);
}
@ -671,26 +578,28 @@ bool Foam::geometricSixDOF::writeData(Ostream& os) const
os.writeKeyword("type") << tab << type() << token::END_STATEMENT << endl;
// Write data
os.writeKeyword("Xrel") << tab << this->Xrel()
os.writeKeyword("Xrel") << tab << Xrel_
<< token::END_STATEMENT << nl;
os.writeKeyword("U") << tab << this->U() << token::END_STATEMENT << nl;
os.writeKeyword("rotationTensor") << tab << this->toRelative()
os.writeKeyword("U") << tab << U_ << token::END_STATEMENT << nl;
os.writeKeyword("rotationTensor") << tab << rotation_
<< token::END_STATEMENT << nl;
os.writeKeyword("omega") << tab << this->omega()
os.writeKeyword("rotationIncrementTensor") << tab << rotIncrement_
<< token::END_STATEMENT << nl;
os.writeKeyword("omega") << tab << 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;
os.writeKeyword("fixedSurge") << tab << fixedSurge_ <<
token::END_STATEMENT << nl;
os.writeKeyword("fixedSway") << tab << fixedSway_ <<
token::END_STATEMENT << nl;
os.writeKeyword("fixedHeave") << tab << fixedHeave_ <<
token::END_STATEMENT << nl;
os.writeKeyword("fixedRoll") << tab << fixedRoll_ <<
token::END_STATEMENT << nl;
os.writeKeyword("fixedPitch") << tab << fixedPitch_ <<
token::END_STATEMENT << nl;
os.writeKeyword("fixedYaw") << tab << fixedYaw_ <<
token::END_STATEMENT << nl << endl;
return os.good();
}

View file

@ -79,24 +79,22 @@ class geometricSixDOF
//- Velocity of mass centroid
dimensionedVector U_;
//- Average velocity of mass centroid at previous time-step
//- Average velocity of mass centroid (evaluated at midstep)
dimensionedVector Uaverage_;
//- Finite rotation. Note: changed during solution process
mutable finiteRotation rotation_;
//- Rotation tensor
tensor rotation_;
//- Rotational increment tensor
tensor rotIncrement_;
//- Rotational velocity about mass centroid
dimensionedVector omega_;
// Average variables that need to be stored
//- Average rotational velocity in relative coordinate system
// (evaluated at midstep)
dimensionedVector omegaAverage_;
//- Average rotational velocity in absolute coordinate system
dimensionedVector omegaAverageAbsolute_;
// ODE controls
@ -107,7 +105,7 @@ class geometricSixDOF
scalarField coeffs_;
//- Motion constraints (given as fixed motion components)
// Motion constraints (given as fixed motion components)
//- Fixed surge (x-translation)
Switch fixedSurge_;
@ -236,7 +234,7 @@ public:
virtual const dimensionedVector& omega() const;
// Displacement and rotation in the absolute coordinate system
// Displacement and velocity in the absolute coordinate system
//- Return position of origin in absolute coordinate system
virtual dimensionedVector X() const;
@ -244,56 +242,22 @@ public:
//- Return velocity of origin
virtual const dimensionedVector& U() const;
//- Return average velocity of origin for the previous time-step
//- Return average velocity of origin (evaluated at midstep)
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
// system (evaluated at midstep)
virtual const dimensionedVector& omegaAverage() const;
//- Return average rotational velocity in absolute coordinate
// system
virtual const dimensionedVector& omegaAverageAbsolute() const;
// Rotations
@ -360,6 +324,7 @@ public:
//- Update ODE after the solution, advancing by delta
virtual void update(const scalar delta);
// Write controls
//- WriteData member function required by regIOobject

View file

@ -198,15 +198,15 @@ Foam::quaternionSixDOF::quaternionSixDOF(const IOobject& io)
Xrel_(dict().lookup("Xrel")),
U_(dict().lookup("U")),
Uaverage_(U_),
Uaverage_("Uaverage", U_),
rotation_
(
vector(dict().lookup("rotationVector")),
dimensionedScalar(dict().lookup("rotationAngle")).value()
),
omega_(dict().lookup("omega")),
omegaAverage_(omega_),
omegaAverageAbsolute_(omega_),
omegaAverage_("omegaAverage", omega_),
omegaAverageAbsolute_("omegaAverageAbsolute", omega_),
coeffs_(13, 0.0),
@ -347,89 +347,13 @@ const Foam::dimensionedVector& Foam::quaternionSixDOF::U() const
return U_;
}
const Foam::dimensionedVector& Foam::quaternionSixDOF::Uaverage() const
{
return Uaverage_;
}
const Foam::finiteRotation& Foam::quaternionSixDOF::rotation() const
{
return rotation_;
}
Foam::vector Foam::quaternionSixDOF::rotVector() const
{
return rotation_.rotVector();
}
Foam::dimensionedScalar Foam::quaternionSixDOF::rotAngle() const
{
return dimensionedScalar("rotAngle", dimless, rotation_.rotAngle());
}
void Foam::quaternionSixDOF::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::quaternionSixDOF::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::quaternionSixDOF::setRotation(const HamiltonRodriguezRot& rot)
{
// Cannot call update rotation: simply set up coefficients
// HJ, 23/Mar/2015
coeffs_[9] = rot.e0();
coeffs_[10] = rot.e1();
coeffs_[11] = rot.e2();
coeffs_[12] = rot.e3();
}
void Foam::quaternionSixDOF::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::quaternionSixDOF::setReferentRotation
(
const HamiltonRodriguezRot& rot
)
{
referentRotation_ = rot;
referentMotionConstraints_ = true;
}
void Foam::quaternionSixDOF::setState(const sixDOFODE& sd)
{
// First set the state in base class subobject
@ -462,25 +386,12 @@ void Foam::quaternionSixDOF::setState(const sixDOFODE& sd)
}
Foam::vector Foam::quaternionSixDOF::rotVectorAverage() const
{
return rotation_.rotVectorAverage();
}
const Foam::dimensionedVector& Foam::quaternionSixDOF::omegaAverage() const
{
return omegaAverage_;
}
const Foam::dimensionedVector&
Foam::quaternionSixDOF::omegaAverageAbsolute() const
{
return omegaAverageAbsolute_;
}
Foam::tensor Foam::quaternionSixDOF::toRelative() const
{
return rotation_.eCurrent().R();
@ -642,14 +553,15 @@ bool Foam::quaternionSixDOF::writeData(Ostream& os) const
os.writeKeyword("type") << tab << type() << token::END_STATEMENT << endl;
// Write data
os.writeKeyword("Xrel") << tab << this->Xrel()
os.writeKeyword("Xrel") << tab << Xrel_
<< token::END_STATEMENT << nl;
os.writeKeyword("U") << tab << this->U() << token::END_STATEMENT << nl;
os.writeKeyword("rotationVector") << tab << this->rotVector()
os.writeKeyword("U") << tab << U_ << token::END_STATEMENT << nl;
os.writeKeyword("rotationVector") << tab << rotation_.rotVector()
<< token::END_STATEMENT << nl;
os.writeKeyword("rotationAngle") << tab << this->rotAngle()
os.writeKeyword("rotationAngle") << tab
<< dimensionedScalar("rotAngle", dimless, rotation_.rotAngle())
<< token::END_STATEMENT << nl;
os.writeKeyword("omega") << tab << this->omega()
os.writeKeyword("omega") << tab << omega_
<< token::END_STATEMENT << nl << nl;
os.writeKeyword("fixedSurge") << tab << fixedSurge_ <<

View file

@ -58,7 +58,7 @@ class quaternionSixDOF
{
// Private data
// Initial body state variables
// Body state variables
//- Displacement relative to spring equilibrium
dimensionedVector Xrel_;
@ -66,7 +66,7 @@ class quaternionSixDOF
//- Velocity of mass centroid
dimensionedVector U_;
//- Average velocity of mass centroid at previous time-step
//- Average velocity of mass centroid (evaluated at midstep)
dimensionedVector Uaverage_;
//- Finite rotation
@ -75,10 +75,8 @@ class quaternionSixDOF
//- Rotational velocity about mass centroid
dimensionedVector omega_;
// Average variables that need to be stored
//- Average rotational velocity in relative coordinate system
// (evaluated at midstep)
dimensionedVector omegaAverage_;
//- Average rotational velocity in absolute coordinate system
@ -89,7 +87,7 @@ class quaternionSixDOF
scalarField coeffs_;
//- Motion constraints (given as fixed motion components)
// Motion constraints (given as fixed motion components)
//- Fixed surge (x-translation)
Switch fixedSurge_;
@ -201,7 +199,7 @@ public:
virtual const dimensionedVector& omega() const;
// Displacement and rotation in the absolute coordinate system
// Displacement and velocity in the absolute coordinate system
//- Return position of origin in absolute coordinate system
virtual dimensionedVector X() const;
@ -209,56 +207,22 @@ public:
//- Return velocity of origin
virtual const dimensionedVector& U() const;
//- Return average velocity of origin for the previous time-step
//- Return average velocity of origin (evaluated at midstep)
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
// system (evaluated at midstep)
virtual const dimensionedVector& omegaAverage() const;
//- Return average rotational velocity in absolute coordinate
// system
virtual const dimensionedVector& omegaAverageAbsolute() const;
// Rotations
@ -325,6 +289,7 @@ public:
//- Update ODE after the solution, advancing by delta
virtual void update(const scalar delta);
// Write controls
//- WriteData member function required by regIOobject

View file

@ -48,7 +48,7 @@ defineRunTimeSelectionTable(sixDOFODE, dictionary);
}
// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
void Foam::sixDOFODE::aitkensRelaxation
(

View file

@ -135,6 +135,8 @@ class sixDOFODE
// Private Member Functions
// Copy control
//- Disallow default bitwise copy construct
sixDOFODE(const sixDOFODE&);
@ -142,9 +144,7 @@ class sixDOFODE
void operator=(const sixDOFODE&);
protected:
// Protected Member Functions
// Aitkens relaxation helper function
//- Update Aitkens relaxation parameters
void aitkensRelaxation(const scalar min, const scalar max);
@ -182,6 +182,7 @@ public:
//- Return autoPtr to the selected sixDOFODE
static autoPtr<sixDOFODE> New(const IOobject& io);
// Destructor
virtual ~sixDOFODE();
@ -271,7 +272,7 @@ public:
virtual const dimensionedVector& omega() const = 0;
// Displacement and rotation in the absolute coordinate system
// Displacement and velocity in the absolute coordinate system
//- Return position of origin in absolute coordinate system
virtual dimensionedVector X() const = 0;
@ -279,57 +280,22 @@ public:
//- Return velocity of origin
virtual const dimensionedVector& U() const = 0;
//- Return average velocity of origin for the previous time-step
//- Return average velocity of origin (evaluated at midstep)
virtual const dimensionedVector& Uaverage() const = 0;
//- Return finite rotation
virtual const finiteRotation& rotation() const = 0;
//- Return rotational vector of body
virtual vector rotVector() const = 0;
//- Return rotation angle of body
virtual dimensionedScalar rotAngle() const = 0;
// Non-access control for setting state variables
//- Set position of origin
virtual void setXrel(const vector& x) = 0;
//- Set velocity of origin
virtual void setU(const vector& u) = 0;
//- Set rotational angle in relative coordinate system
virtual void setRotation(const HamiltonRodriguezRot& rot) = 0;
//- Set rotational velocity in relative coordinate system
virtual void setOmega(const vector& omega) = 0;
//- Set referent coordinate system to apply constraints
virtual void setReferentRotation
(
const HamiltonRodriguezRot& rot
) = 0;
//- 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 = 0;
//- Return average rotational velocity in relative coordinate
// system
// system (evaluated at midstep)
virtual const dimensionedVector& omegaAverage() const = 0;
//- Return average rotational velocity in absolute coordinate
// system
virtual const dimensionedVector&
omegaAverageAbsolute() const = 0;
// Rotations