ENH: Using proper expression for accelerating frame terms.
No "axis" required - when considering omegaDot the concept becomes meaningless. Centre of rotation can be specified, to avoid it needing to be the origin, and to allow it to move.
This commit is contained in:
graham
parent
5528940ae5
commit
cbc4cee935
@ -22,24 +22,24 @@
|
||||
}
|
||||
|
||||
|
||||
IOobject angularVelocityMagHeader
|
||||
IOobject angularVelocityHeader
|
||||
(
|
||||
"angularVelocityMag",
|
||||
"angularVelocity",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
);
|
||||
|
||||
autoPtr<uniformDimensionedScalarField> angularVelocityMagPtr;
|
||||
autoPtr<uniformDimensionedVectorField> angularVelocityPtr;
|
||||
|
||||
if (angularVelocityMagHeader.headerOk())
|
||||
if (angularVelocityHeader.headerOk())
|
||||
{
|
||||
Info<< " Reading " << angularVelocityMagHeader.name() << endl;
|
||||
Info<< " Reading " << angularVelocityHeader.name() << endl;
|
||||
|
||||
angularVelocityMagPtr.reset
|
||||
angularVelocityPtr.reset
|
||||
(
|
||||
new uniformDimensionedScalarField(angularVelocityMagHeader)
|
||||
new uniformDimensionedVectorField(angularVelocityHeader)
|
||||
);
|
||||
}
|
||||
|
||||
@ -66,45 +66,23 @@
|
||||
}
|
||||
|
||||
|
||||
IOobject axisHeader
|
||||
IOobject centreOfRotationHeader
|
||||
(
|
||||
"axis",
|
||||
"centreOfRotation",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
);
|
||||
|
||||
autoPtr<uniformDimensionedVectorField> axisPtr;
|
||||
autoPtr<uniformDimensionedVectorField> centreOfRotationPtr;
|
||||
|
||||
if (axisHeader.headerOk())
|
||||
if (centreOfRotationHeader.headerOk())
|
||||
{
|
||||
Info<< " Reading " << axisHeader.name() << endl;
|
||||
Info<< " Reading " << centreOfRotationHeader.name() << endl;
|
||||
|
||||
axisPtr.reset
|
||||
centreOfRotationPtr.reset
|
||||
(
|
||||
new uniformDimensionedVectorField(axisHeader)
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
IOobject axisRefPointHeader
|
||||
(
|
||||
"axisRefPoint",
|
||||
runTime.timeName(),
|
||||
mesh,
|
||||
IOobject::MUST_READ,
|
||||
IOobject::AUTO_WRITE
|
||||
);
|
||||
|
||||
autoPtr<uniformDimensionedVectorField> axisRefPointPtr;
|
||||
|
||||
if (axisRefPointHeader.headerOk())
|
||||
{
|
||||
Info<< " Reading " << axisRefPointHeader.name() << endl;
|
||||
|
||||
axisRefPointPtr.reset
|
||||
(
|
||||
new uniformDimensionedVectorField(axisRefPointHeader)
|
||||
new uniformDimensionedVectorField(centreOfRotationHeader)
|
||||
);
|
||||
}
|
||||
|
||||
@ -47,15 +47,15 @@ Foam::NonInertialFrameForce<CloudType>::NonInertialFrameForce
|
||||
)
|
||||
),
|
||||
W_(vector::zero),
|
||||
omegaMagName_
|
||||
omegaName_
|
||||
(
|
||||
this->coeffs().template lookupOrDefault<word>
|
||||
(
|
||||
"angularVelocityMagName",
|
||||
"angularVelocityMag"
|
||||
"angularVelocityName",
|
||||
"angularVelocity"
|
||||
)
|
||||
),
|
||||
omegaMag_(0.0),
|
||||
omega_(vector::zero),
|
||||
omegaDotName_
|
||||
(
|
||||
this->coeffs().template lookupOrDefault<word>
|
||||
@ -65,25 +65,15 @@ Foam::NonInertialFrameForce<CloudType>::NonInertialFrameForce
|
||||
)
|
||||
),
|
||||
omegaDot_(vector::zero),
|
||||
axisName_
|
||||
centreOfRotationName_
|
||||
(
|
||||
this->coeffs().template lookupOrDefault<word>
|
||||
(
|
||||
"axisName",
|
||||
"axis"
|
||||
"centreOfRotationName",
|
||||
"centreOfRotation"
|
||||
)
|
||||
),
|
||||
axis_(vector::zero),
|
||||
hasAxis_(false),
|
||||
axisRefPointName_
|
||||
(
|
||||
this->coeffs().template lookupOrDefault<word>
|
||||
(
|
||||
"axisRefPointName",
|
||||
"axisRefPoint"
|
||||
)
|
||||
),
|
||||
axisRefPoint_(vector::zero)
|
||||
centreOfRotation_(vector::zero)
|
||||
{}
|
||||
|
||||
|
||||
@ -96,15 +86,12 @@ Foam::NonInertialFrameForce<CloudType>::NonInertialFrameForce
|
||||
ParticleForce<CloudType>(niff),
|
||||
WName_(niff.WName_),
|
||||
W_(niff.W_),
|
||||
omegaMagName_(niff.omegaMagName_),
|
||||
omegaMag_(niff.omegaMag_),
|
||||
omegaName_(niff.omegaName_),
|
||||
omega_(niff.omega_),
|
||||
omegaDotName_(niff.omegaDotName_),
|
||||
omegaDot_(niff.omegaDot_),
|
||||
axisName_(niff.axisName_),
|
||||
axis_(niff.axis_),
|
||||
hasAxis_(niff.hasAxis_),
|
||||
axisRefPointName_(niff.axisRefPointName_),
|
||||
axisRefPoint_(niff.axisRefPoint_)
|
||||
centreOfRotationName_(niff.centreOfRotationName_),
|
||||
centreOfRotation_(niff.centreOfRotation_)
|
||||
{}
|
||||
|
||||
|
||||
@ -121,71 +108,18 @@ template<class CloudType>
|
||||
void Foam::NonInertialFrameForce<CloudType>::cacheFields(const bool store)
|
||||
{
|
||||
W_ = vector::zero;
|
||||
omegaMag_ = 0.0;
|
||||
omega_ = vector::zero;
|
||||
omegaDot_ = vector::zero;
|
||||
axis_ = vector::zero;
|
||||
hasAxis_ = false;
|
||||
axisRefPoint_ = vector::zero;
|
||||
centreOfRotation_ = vector::zero;
|
||||
|
||||
if (store)
|
||||
{
|
||||
if
|
||||
(
|
||||
this->mesh().template
|
||||
foundObject<uniformDimensionedScalarField>(omegaMagName_)
|
||||
)
|
||||
{
|
||||
uniformDimensionedScalarField omegaMag = this->mesh().template
|
||||
lookupObject<uniformDimensionedScalarField>(omegaMagName_);
|
||||
|
||||
omegaMag_ = omegaMag.value();
|
||||
|
||||
// If omegaMag is found, require that the axis and axisRefPoint is
|
||||
// found.
|
||||
uniformDimensionedVectorField a = this->mesh().template
|
||||
lookupObject<uniformDimensionedVectorField>(axisName_);
|
||||
|
||||
axis_ = a.value();
|
||||
|
||||
hasAxis_ = true;
|
||||
|
||||
scalar axisMag = mag(axis_);
|
||||
|
||||
if (mag(axis_) < SMALL)
|
||||
{
|
||||
FatalErrorIn
|
||||
(
|
||||
"void Foam::NonInertialFrameForce<CloudType>::"
|
||||
"cacheFields(const bool store)"
|
||||
) << axisName_ << " " << axis_ << " too small."
|
||||
<< abort(FatalError);
|
||||
}
|
||||
|
||||
axis_ /= axisMag;
|
||||
|
||||
uniformDimensionedVectorField axisRefPoint = this->mesh().template
|
||||
lookupObject<uniformDimensionedVectorField>(axisRefPointName_);
|
||||
|
||||
axisRefPoint_ = axisRefPoint.value();
|
||||
|
||||
// Only look for omegaDot is omegaMag is found, optional.
|
||||
if
|
||||
this->mesh().template foundObject<uniformDimensionedVectorField>
|
||||
(
|
||||
this->mesh().template
|
||||
foundObject<uniformDimensionedVectorField>(omegaDotName_)
|
||||
WName_
|
||||
)
|
||||
{
|
||||
uniformDimensionedVectorField omegaDot = this->mesh().template
|
||||
lookupObject<uniformDimensionedVectorField>(omegaDotName_);
|
||||
|
||||
omegaDot_ = omegaDot.value();
|
||||
}
|
||||
}
|
||||
|
||||
if
|
||||
(
|
||||
this->mesh().template
|
||||
foundObject<uniformDimensionedVectorField>(WName_)
|
||||
)
|
||||
{
|
||||
uniformDimensionedVectorField W = this->mesh().template
|
||||
@ -193,13 +127,50 @@ void Foam::NonInertialFrameForce<CloudType>::cacheFields(const bool store)
|
||||
|
||||
W_ = W.value();
|
||||
}
|
||||
else if (!hasAxis_)
|
||||
{
|
||||
WarningIn
|
||||
|
||||
if
|
||||
(
|
||||
this->mesh().template foundObject<uniformDimensionedVectorField>
|
||||
(
|
||||
"void Foam::NonInertialFrameForce<CloudType>::"
|
||||
"cacheFields(const bool store)"
|
||||
) << "No " << typeName << " variables found." << endl;
|
||||
omegaName_
|
||||
)
|
||||
)
|
||||
{
|
||||
uniformDimensionedVectorField omega = this->mesh().template
|
||||
lookupObject<uniformDimensionedVectorField>(omegaName_);
|
||||
|
||||
omega_ = omega.value();
|
||||
}
|
||||
|
||||
if
|
||||
(
|
||||
this->mesh().template foundObject<uniformDimensionedVectorField>
|
||||
(
|
||||
omegaDotName_
|
||||
)
|
||||
)
|
||||
{
|
||||
uniformDimensionedVectorField omegaDot = this->mesh().template
|
||||
lookupObject<uniformDimensionedVectorField>(omegaDotName_);
|
||||
|
||||
omegaDot_ = omegaDot.value();
|
||||
}
|
||||
|
||||
if
|
||||
(
|
||||
this->mesh().template foundObject<uniformDimensionedVectorField>
|
||||
(
|
||||
centreOfRotationName_
|
||||
)
|
||||
)
|
||||
{
|
||||
uniformDimensionedVectorField omegaDot = this->mesh().template
|
||||
lookupObject<uniformDimensionedVectorField>
|
||||
(
|
||||
centreOfRotationName_
|
||||
);
|
||||
|
||||
centreOfRotation_ = omegaDot.value();
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -217,24 +188,16 @@ Foam::forceSuSp Foam::NonInertialFrameForce<CloudType>::calcNonCoupled
|
||||
{
|
||||
forceSuSp value(vector::zero, 0.0);
|
||||
|
||||
value.Su() += -mass*W_;
|
||||
const vector& r = p.position() - centreOfRotation_;
|
||||
|
||||
if (hasAxis_)
|
||||
{
|
||||
const vector pRel = p.position() - axisRefPoint_;
|
||||
|
||||
const vector r = pRel - axis_*(axis_ & pRel);
|
||||
|
||||
vector omega = axis_*omegaMag_;
|
||||
|
||||
value.Su() +=
|
||||
mass
|
||||
*(
|
||||
(r ^ omegaDot_)
|
||||
+ 2.0*(p.U() ^ omega)
|
||||
+ (omega ^ (r ^ omega))
|
||||
);
|
||||
}
|
||||
value.Su() =
|
||||
mass
|
||||
*(
|
||||
-W_
|
||||
+ (r ^ omegaDot_)
|
||||
+ 2.0*(p.U() ^ omega_)
|
||||
+ (omega_ ^ (r ^ omega_))
|
||||
);
|
||||
|
||||
return value;
|
||||
}
|
||||
|
||||
@ -62,36 +62,23 @@ class NonInertialFrameForce
|
||||
//- The linear acceleration of the reference frame
|
||||
vector W_;
|
||||
|
||||
//- Name of the angular velocity magnitude field
|
||||
word omegaMagName_;
|
||||
//- Name of the angular velocity field
|
||||
word omegaName_;
|
||||
|
||||
//- The magnitude of the angular velocity of the reference frame,
|
||||
// combines with axis to give omega
|
||||
scalar omegaMag_;
|
||||
//- The angular velocity of the reference frame
|
||||
vector omega_;
|
||||
|
||||
//- Name of the angular acceleration field
|
||||
word omegaDotName_;
|
||||
|
||||
//- Pointer to the angular acceleration of the reference frame
|
||||
//- The angular acceleration of the reference frame
|
||||
vector omegaDot_;
|
||||
|
||||
//- Name of the axis field
|
||||
word axisName_;
|
||||
//- Name of the centre of rotation field
|
||||
word centreOfRotationName_;
|
||||
|
||||
//- Pointer to the axis of rotation - assumed to be a unit vector.
|
||||
// duplication of omega to allow situations where omega = 0 and
|
||||
// omegaDot != 0.
|
||||
vector axis_;
|
||||
|
||||
// Boolean flag for whether rotational motion is active
|
||||
bool hasAxis_;
|
||||
|
||||
//- Name of the axisRefPoint field
|
||||
word axisRefPointName_;
|
||||
|
||||
//- Pointer to a point in non-inertial space on the axis of rotation
|
||||
// (omega), used to calculate r.
|
||||
vector axisRefPoint_;
|
||||
//- The centre of rotation of the reference frame
|
||||
vector centreOfRotation_;
|
||||
|
||||
|
||||
public:
|
||||
@ -141,14 +128,9 @@ public:
|
||||
//- Return the angular acceleration of the reference frame
|
||||
inline const vector& omegaDot() const;
|
||||
|
||||
//- Return the axis of rotation
|
||||
inline const vector& axis() const;
|
||||
//- Return the centre of rotation of the reference frame
|
||||
inline const vector& centreOfRotation() const;
|
||||
|
||||
//- Return bool stating if the frame is rotating
|
||||
inline bool hasAxis() const;
|
||||
|
||||
//- Return the point in non-inertial space on the axis of rotation
|
||||
inline const vector& axisRefPoint() const;
|
||||
|
||||
// Evaluation
|
||||
|
||||
|
||||
@ -35,7 +35,7 @@ inline const Foam::vector& Foam::NonInertialFrameForce<CloudType>::W() const
|
||||
template<class CloudType>
|
||||
inline const Foam::vector& Foam::NonInertialFrameForce<CloudType>::omega() const
|
||||
{
|
||||
return omegaMag_*axis_;
|
||||
return omega_;
|
||||
}
|
||||
|
||||
|
||||
@ -49,24 +49,9 @@ Foam::NonInertialFrameForce<CloudType>::omegaDot() const
|
||||
|
||||
template<class CloudType>
|
||||
inline const Foam::vector&
|
||||
Foam::NonInertialFrameForce<CloudType>::axis() const
|
||||
Foam::NonInertialFrameForce<CloudType>::centreOfRotation() const
|
||||
{
|
||||
return axis_;
|
||||
}
|
||||
|
||||
|
||||
template<class CloudType>
|
||||
inline bool Foam::NonInertialFrameForce<CloudType>::hasAxis() const
|
||||
{
|
||||
return hasAxis_;
|
||||
}
|
||||
|
||||
|
||||
template<class CloudType>
|
||||
inline const Foam::vector&
|
||||
Foam::NonInertialFrameForce<CloudType>::axisRefPoint() const
|
||||
{
|
||||
return axisRefPoint_;
|
||||
return centreOfRotation_;
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -92,9 +92,8 @@ Foam::forceSuSp Foam::SRFForce<CloudType>::calcNonCoupled
|
||||
const typename SRF::SRFModel& srf = *srfPtr_;
|
||||
|
||||
const vector& omega = srf.omega().value();
|
||||
const vector& axis = srf.axis();
|
||||
|
||||
const vector r = p.position() - axis*(axis & p.position());
|
||||
const vector& r = p.position();
|
||||
|
||||
// Coriolis and centrifugal acceleration terms
|
||||
value.Su() = mass*(2.0*(p.U() ^ omega) + (omega ^ (r ^ omega)));
|
||||
|
||||
Loading…
Reference in New Issue
Block a user