Updates for multiple solution of 6-DOF ODE

This commit is contained in:
Hrvoje Jasak 2015-04-07 13:29:39 +01:00 committed by Dominik Christ
parent 6f424c2635
commit 736f1b496a
7 changed files with 157 additions and 48 deletions

View file

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

View file

@ -52,7 +52,7 @@ class finiteRotation
// Private data
//- Initial rotation
const HamiltonRodriguezRot eInitial_;
HamiltonRodriguezRot eInitial_;
//- Rotational tensor
tensor rotTensor_;
@ -63,13 +63,6 @@ class finiteRotation
// Private Member Functions
//- Disallow default bitwise copy construct
finiteRotation(const finiteRotation&);
//- Disallow default bitwise assignment
void operator=(const finiteRotation&);
//- Calculate unit rotation vector from given rotation tensor
static vector rotVector(const tensor& rotT);

View file

@ -74,17 +74,17 @@ class HamiltonRodriguezRot
//- Calculate R_ - inertial to body coord. sys. rotation
inline void calculateR() const
{
R_.xx() = 2 * (e1_*e1_ + e0_*e0_ - 0.5);
R_.xy() = 2 * (e1_*e2_ + e0_*e3_);
R_.xz() = 2 * (e1_*e3_ - e0_*e2_);
R_.xx() = 2*(e1_*e1_ + e0_*e0_ - 0.5);
R_.xy() = 2*(e1_*e2_ + e0_*e3_);
R_.xz() = 2*(e1_*e3_ - e0_*e2_);
R_.yx() = 2 * (e2_*e1_ - e0_*e3_);
R_.yy() = 2 * (e2_*e2_ + e0_*e0_ - 0.5);
R_.yz() = 2 * (e2_*e3_ + e0_*e1_);
R_.yx() = 2*(e2_*e1_ - e0_*e3_);
R_.yy() = 2*(e2_*e2_ + e0_*e0_ - 0.5);
R_.yz() = 2*(e2_*e3_ + e0_*e1_);
R_.zx() = 2 * (e3_*e1_ + e0_*e2_);
R_.zy() = 2 * (e3_*e2_ - e0_*e1_);
R_.zz() = 2 * (e3_*e3_ + e0_*e0_ - 0.5);
R_.zx() = 2*(e3_*e1_ + e0_*e2_);
R_.zy() = 2*(e3_*e2_ - e0_*e1_);
R_.zz() = 2*(e3_*e3_ + e0_*e0_ - 0.5);
}
//- Calculate Gt()

View file

@ -30,7 +30,7 @@ Description
Author
Dubravko Matijasevic, FSB Zagreb. All rights reserved.
\*----------------------------------------------------------------------------*/
\*---------------------------------------------------------------------------*/
#include "objectRegistry.H"
#include "sixDOFbodies.H"

View file

@ -35,17 +35,6 @@ Author
#include "objectRegistry.H"
#include "sixDOFqODE.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
// Runtime type information
// Not possible because of I/O error: incorrect type, expecting dictionary
// HJ, 11/Feb/2008
// namespace Foam
// {
// defineTypeNameAndDebug(sixDOFqODE, 0);
// }
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::sixDOFqODE::setCoeffs()
@ -155,6 +144,49 @@ Foam::sixDOFqODE::sixDOFqODE(const IOobject& io)
}
// Construct as copy
Foam::sixDOFqODE::sixDOFqODE
(
const word& name,
const sixDOFqODE& sd
)
:
IOdictionary
(
IOobject
(
name,
sd.instance(),
sd.local(),
sd.db(),
IOobject::NO_READ,
IOobject::NO_WRITE
)
),
mass_(sd.mass_.name(), sd.mass_),
momentOfInertia_(sd.momentOfInertia_.name(), sd.momentOfInertia_),
Xequilibrium_(sd.Xequilibrium_.name(), sd.Xequilibrium_),
linSpringCoeffs_(sd.linSpringCoeffs_.name(), sd.linSpringCoeffs_),
linDampingCoeffs_(sd.linDampingCoeffs_.name(), sd.linDampingCoeffs_),
Xrel_(sd.Xrel_.name(), sd.Xrel_),
U_(sd.U_.name(), sd.U_),
Uaverage_(sd.Uaverage_.name(), sd.Uaverage_),
rotation_(sd.rotation_),
omega_(sd.omega_.name(), sd.omega_),
omegaAverage_(sd.omegaAverage_.name(), sd.omegaAverage_),
omegaAverageAbsolute_
(
sd.omegaAverageAbsolute_.name(),
sd.omegaAverageAbsolute_
),
force_(sd.force_.name(), sd.force_),
moment_(sd.moment_.name(), sd.moment_),
forceRelative_(sd.forceRelative_.name(), sd.forceRelative_),
momentRelative_(sd.momentRelative_.name(), sd.momentRelative_),
coeffs_(sd.coeffs_)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::sixDOFqODE::~sixDOFqODE()
@ -163,6 +195,31 @@ Foam::sixDOFqODE::~sixDOFqODE()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::sixDOFqODE::setState(const sixDOFqODE& sd)
{
mass_ = sd.mass_;
momentOfInertia_ = sd.momentOfInertia_;
Xequilibrium_ = sd.Xequilibrium_;
linSpringCoeffs_ = sd.linSpringCoeffs_;
linDampingCoeffs_ = sd.linDampingCoeffs_;
Xrel_ = sd.Xrel_;
U_ = sd.U_;
Uaverage_ = sd.Uaverage_;
rotation_ = sd.rotation_;
omega_ = sd.omega_;
omegaAverage_ = sd.omegaAverage_;
omegaAverageAbsolute_ = sd.omegaAverageAbsolute_;
force_ = sd.force_;
moment_ = sd.moment_;
forceRelative_ = sd.forceRelative_;
momentRelative_ = sd.momentRelative_;
// Copy ODE coefficients: this carries actual ODE state
// HJ, 23/Mar/2015
coeffs_ = sd.coeffs_;
}
void Foam::sixDOFqODE::derivatives
(
const scalar x,

View file

@ -173,18 +173,18 @@ class sixDOFqODE
public:
// //- Runtime type information
// Not possible because of I/O error: incorrect type, expecting dictionary
// HJ, 11/Feb/2008
// TypeName("sixDOFqODE");
// Constructors
//- Construct from dictionary
sixDOFqODE(const IOobject& io);
//- Construct sixDOFqODE, changing name
sixDOFqODE
(
const word& name,
const sixDOFqODE& sd
);
// Destructor
@ -246,14 +246,21 @@ public:
inline dimensionedScalar rotAngle() const;
// Non-constant access to velocities
// Non-constant access to
//- Return access to velocity of origin
inline dimensionedVector& U();
//- Set position of origin
inline void setXrel(const vector& x);
//- Return access to rotational velocity in relative
//- Set velocity of origin
inline void setU(const vector& u);
//- Set rotational angle in relative
// coordinate system
inline dimensionedVector& omega();
inline void setRotation(const HamiltonRodriguezRot& rot);
//- Set rotational velocity in relative
// coordinate system
inline void setOmega(const vector& omega);
// Average motion per time-step
@ -261,10 +268,12 @@ public:
//- Return average rotational vector of body
inline vector rotVectorAverage() const;
//- Return average rotational velocity in relative coordinate system
//- Return average rotational velocity in
// relative coordinate system
inline const dimensionedVector& omegaAverage() const;
//- Return average rotational velocity in absolute coordinate system
//- Return average rotational velocity in
// absolute coordinate system
inline const dimensionedVector& omegaAverageAbsolute() const;
@ -313,6 +322,9 @@ public:
//- Return transformation tensor between new and previous rotation
inline const tensor& rotIncrementTensor() const;
//- Set ODE parameters from another ODE
void setState(const sixDOFqODE&);
// ODE parameters

View file

@ -118,15 +118,51 @@ Foam::dimensionedScalar Foam::sixDOFqODE::rotAngle() const
}
Foam::dimensionedVector& Foam::sixDOFqODE::U()
void Foam::sixDOFqODE::setXrel(const vector& x)
{
return U_;
Xrel_.value() = x;
// Set X in coefficients
coeffs_[0] = x.x();
coeffs_[1] = x.y();
coeffs_[2] = x.z();
}
Foam::dimensionedVector& Foam::sixDOFqODE::omega()
void Foam::sixDOFqODE::setU(const vector& U)
{
return omega_;
U_.value() = U;
Uaverage_ = U_;
// Set U in coefficients
coeffs_[3] = U.x();
coeffs_[4] = U.y();
coeffs_[5] = U.z();
}
void Foam::sixDOFqODE::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::sixDOFqODE::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();
}