songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'feature-sixdof-constraints-restraints' into 'develop'

Browse Source 
ENH: sixDoF: add max/min angle constraints to the axis constraint

See merge request Development/openfoam!613
This commit is contained in:
Andrew Heather 2023-06-14 14:30:55 +00:00
commit d2e8c75dc9
12 changed files with 401 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

257
src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/axis/sixDoFRigidBodyMotionAxisConstraint.C
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2015 OpenFOAM Foundation
Copyright (C) 2023 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -28,6 +29,7 @@ License
#include "sixDoFRigidBodyMotionAxisConstraint.H"
#include "addToRunTimeSelectionTable.H"
#include "sixDoFRigidBodyMotion.H"
#include "unitConversion.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -47,6 +49,42 @@ namespace sixDoFRigidBodyMotionConstraints
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::label Foam::sixDoFRigidBodyMotionConstraints::axis::rotationSector
(
const vector& oldDir,
const vector& newDir
) const
{
const scalar thetaDir = (oldDir ^ newDir) & axis_;
if (equal(thetaDir, 0))
{
return 0;
}
return label(sign(thetaDir));
}
bool Foam::sixDoFRigidBodyMotionConstraints::axis::calcDir
(
const vector& fm,
const bool rotationSector
) const
{
const scalar fmDir = axis_ & fm;
if (equal(fmDir, 0))
{
return rotationSector;
}
return (label(sign(fmDir)) == 1) ? true : false;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionConstraints::axis::axis
@ -57,18 +95,16 @@ Foam::sixDoFRigidBodyMotionConstraints::axis::axis
)
:
sixDoFRigidBodyMotionConstraint(name, sDoFRBMCDict, motion),
axis_()
refQ_(),
axis_(),
maxCWThetaPtr_(),
maxCCWThetaPtr_(),
degrees_(false)
{
read(sDoFRBMCDict);
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotionConstraints::axis::~axis()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::sixDoFRigidBodyMotionConstraints::axis::constrainTranslation
@ -83,7 +119,115 @@ void Foam::sixDoFRigidBodyMotionConstraints::axis::constrainRotation
pointConstraint& pc
) const
{
pc.combine(pointConstraint(Tuple2<label, vector>(2, axis_)));
if (!(maxCWThetaPtr_ && maxCCWThetaPtr_))
{
pc.combine(pointConstraint(Tuple2<label, vector>(2, axis_)));
return;
}
// Calculate principal directions of the body
const vector refDir
(
rotationTensor(vector(1, 0 ,0), axis_) & vector(0, 1, 0)
);
const vector oldDir
(
(refQ_ & refDir).removeCollinear(axis_).normalise()
);
const vector newDir
(
(motion().orientation() & refDir).removeCollinear(axis_).normalise()
);
// Find the index of the rotation sector that the body resides
const label rotationSectorIndex = rotationSector(oldDir, newDir);
if (!rotationSectorIndex)
{
// The body resides at the reference orientation
pc.combine(pointConstraint(Tuple2<label, vector>(2, axis_)));
return;
}
const bool rotationSector = (rotationSectorIndex == 1) ? true : false;
// Calculate the directions of total momentum and force acting on the body
const bool angularMomentumDir =
calcDir
(
motion().state().pi(),
rotationSector
);
const bool torqueDir = calcDir(motion().state().tau(), rotationSector);
// Calculate the rotation angle of the body wrt the reference orientation
const scalar theta = mag(acos(min(oldDir & newDir, scalar(1))));
// Calculate maximum clockwise and counterclockwise rotation angles
const scalar t = motion().time().timeOutputValue();
const scalar maxCWTheta =
degrees_
? mag(degToRad(maxCWThetaPtr_->value(t)))
: mag(maxCWThetaPtr_->value(t));
const scalar maxCCWTheta =
degrees_
? mag(degToRad(maxCCWThetaPtr_->value(t)))
: mag(maxCCWThetaPtr_->value(t));
// Apply the constraints according to various conditions
if
(
(rotationSector && (theta < maxCCWTheta))
|| (!rotationSector && (theta < maxCWTheta))
)
{
pc.combine(pointConstraint(Tuple2<label, vector>(2, axis_)));
}
else
{
if (rotationSector == angularMomentumDir)
{
const scalar magPi = mag(motion().state().pi());
if (equal(magPi, scalar(0)) && rotationSector != torqueDir)
{
pc.combine(pointConstraint(Tuple2<label, vector>(2, axis_)));
}
else
{
// Constrain all rotational motions
pc.combine(pointConstraint(Tuple2<label, vector>(3, Zero)));
}
}
else
{
// If there is a difference between the directions of
// body rotation and of torque, release the body
pc.combine(pointConstraint(Tuple2<label, vector>(2, axis_)));
}
}
if (motion().report())
{
Info
<< " old direction = " << oldDir << nl
<< " new direction = " << newDir << nl
<< " rotationSector = " << rotationSector << nl
<< " theta = " << sign((oldDir ^ newDir) & axis_)*theta << nl
<< " torque = " << motion().state().tau() << nl
<< " torque dir = " << torqueDir << nl
<< " angular momentum = " << motion().state().pi() << nl
<< " angular momentum dir = " << angularMomentumDir << nl
<< endl;
}
}
@ -92,21 +236,76 @@ bool Foam::sixDoFRigidBodyMotionConstraints::axis::read
const dictionary& sDoFRBMCDict
)
{
sixDoFRigidBodyMotionConstraint::read(sDoFRBMCDict);
if (!sixDoFRigidBodyMotionConstraint::read(sDoFRBMCDict))
{
return false;
}
sDoFRBMCCoeffs_.readEntry("axis", axis_);
scalar magFixedAxis(mag(axis_));
axis_.normalise();
if (magFixedAxis > VSMALL)
if (mag(axis_) < VSMALL)
{
axis_ /= magFixedAxis;
FatalIOErrorInFunction(sDoFRBMCDict)
<< "The entry 'axis' cannot have zero length: " << axis_
<< exit(FatalIOError);
}
else
if (sDoFRBMCCoeffs_.found("thetaUnits"))
{
FatalErrorInFunction
<< "axis has zero length"
<< abort(FatalError);
const word thetaUnits(sDoFRBMCCoeffs_.getWord("thetaUnits"));
if (thetaUnits == "degrees")
{
degrees_ = true;
}
else if (thetaUnits == "radians")
{
degrees_ = false;
}
else
{
FatalIOErrorInFunction(sDoFRBMCCoeffs_)
<< "The units of thetaUnits can be either degrees or radians"
<< exit(FatalIOError);
}
maxCWThetaPtr_.reset
(
Function1<scalar>::New
(
"maxClockwiseTheta",
sDoFRBMCCoeffs_,
&motion().time()
)
);
maxCCWThetaPtr_.reset
(
Function1<scalar>::New
(
"maxCounterclockwiseTheta",
sDoFRBMCCoeffs_,
&motion().time()
)
);
refQ_ =
sDoFRBMCCoeffs_.getOrDefault<tensor>("referenceOrientation", I);
if (mag(mag(refQ_) - sqrt(3.0)) > ROOTSMALL)
{
FatalIOErrorInFunction(sDoFRBMCCoeffs_)
<< "The entry 'referenceOrientation' " << refQ_
<< " is not a rotation tensor. "
<< "mag(referenceOrientation) - sqrt(3) = "
<< mag(refQ_) - sqrt(3.0) << nl
<< exit(FatalIOError);
}
}
return true;
@ -119,6 +318,32 @@ void Foam::sixDoFRigidBodyMotionConstraints::axis::write
) const
{
os.writeEntry("axis", axis_);
if (maxCWThetaPtr_ && maxCCWThetaPtr_)
{
if (degrees_)
{
os.writeEntry("thetaUnits", "degrees");
}
else
{
os.writeEntry("thetaUnits", "radians");
}
if (maxCWThetaPtr_)
{
maxCWThetaPtr_->writeData(os);
}
if (maxCCWThetaPtr_)
{
maxCCWThetaPtr_->writeData(os);
}
os.writeEntry("referenceOrientation", refQ_);
}
}
// ************************************************************************* //

90
src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/axis/sixDoFRigidBodyMotionAxisConstraint.H
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2014 OpenFOAM Foundation
Copyright (C) 2023 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -27,18 +28,67 @@ Class
Foam::sixDoFRigidBodyMotionConstraints::axis
Description
Orientation constraint:
may only rotate around a fixed axis.
This constraint imposes an orientation limitation where bodies
are restricted to rotate only around a fixed axis, optionally with
the inclusion of maximum and minimum rotation angle constraints.
Usage
Minimal example by \c constant/dynamicMeshDict.sixDoFRigidBodyMotionCoeffs:
\verbatim
constraints
{
constrainRotation1
{
// Mandatory entries
sixDoFRigidBodyMotionConstraint axis;
axis <vector>;
// Optional entries
maxClockwiseTheta <Function1<scalar>>;
maxCounterclockwiseTheta <Function1<scalar>>;
thetaUnits <word>;
referenceOrientation <tensor>;
}
}
\endverbatim
where the entries mean:
\table
Property | Description | Type | Reqd | Deflt
sixDoFRigidBodyMotionConstraint | Type name: axis | word | yes | -
axis | Reference rotation axis fixed in global space <!--
--> | vector | yes | -
maxClockwiseTheta | Maximum clockwise rotation angle about the <!--
--> rotation axis | Function1\<scalar\> | no | -
maxCounterlockwiseTheta | Maximum counterclockwise rotation <!--
--> angle about the rotation axis | Function1\<scalar\> | no | -
thetaUnits | Units of theta angles | word | no | radians
referenceOrientation | Reference orientation where there is no moment <!--
--> | tensor | no | I
\endtable
The inherited entries are elaborated in:
- \link Function1.H \endlink
Notes
- The units for \c thetaUnits can be specified
as either \c degrees or \c radians.
- The \c maxClockwiseTheta and \c maxCounterlockwiseTheta
are always non-negative.
- Negative and positive \c theta correspond to clockwise
and counterclockwise rotation sectors with respect to
the reference orientation, respectively.
SourceFiles
sixDoFRigidBodyMotionAxisConstraint.C
\*---------------------------------------------------------------------------*/
#ifndef sixDoFRigidBodyMotionAxisConstraint_H
#define sixDoFRigidBodyMotionAxisConstraint_H
#ifndef Foam_sixDoFRigidBodyMotionAxisConstraint_H
#define Foam_sixDoFRigidBodyMotionAxisConstraint_H
#include "sixDoFRigidBodyMotionConstraint.H"
#include "Function1.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -56,12 +106,35 @@ class axis
:
public sixDoFRigidBodyMotionConstraint
{
// Private Data
// Private data
//- Reference orientation where there is no moment
tensor refQ_;
//- Reference axis in global space
//- Reference rotation axis fixed in global space
vector axis_;
//- Maximum clockwise rotation angle around the rotation axis
autoPtr<Function1<scalar>> maxCWThetaPtr_;
//- Maximum counterclockwise rotation angle around the rotation axis
autoPtr<Function1<scalar>> maxCCWThetaPtr_;
//- Flag to enable to use maxTheta in degrees rather than in radians
bool degrees_;
// Private Member Functions
//- Return the index of the rotation sector that the
//- body resides with respect to the reference orientation
// -1: clockwise, 0: reference orientation, 1: counterclockwise
label rotationSector(const vector& oldDir, const vector& newDir) const;
//- Return the direction of the given force or momentum
// false: clockwise, true: counterclockwise
bool calcDir(const vector& fm, const bool rotationSector) const;
public:
@ -88,9 +161,12 @@ public:
);
}
//- No copy assignment
void operator=(const sixDoFRigidBodyMotionConstraint&) = delete;
//- Destructor
virtual ~axis();
virtual ~axis() = default;
// Member Functions

7
src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/constraints/sixDoFRigidBodyMotionConstraint/sixDoFRigidBodyMotionConstraint.H
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2014 OpenFOAM Foundation
Copyright (C) 2023 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -162,6 +163,12 @@ public:
return sDoFRBMCCoeffs_;
}
//- Return const access to motion
const sixDoFRigidBodyMotion& motion() const noexcept { return motion_; }
// I-O
//- Write
virtual void write(Ostream&) const;
};

30
src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
Copyright (C) 2016-2020 OpenCFD Ltd.
Copyright (C) 2016-2023 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -90,6 +90,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion(const Time& time)
aRelax_(1.0),
aDamp_(1.0),
report_(false),
updateConstraints_(false),
solver_(nullptr)
{}
@ -130,6 +131,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
aRelax_(dict.getOrDefault<scalar>("accelerationRelaxation", 1)),
aDamp_(dict.getOrDefault<scalar>("accelerationDamping", 1)),
report_(dict.getOrDefault("report", false)),
updateConstraints_(dict.getOrDefault("updateConstraints", false)),
solver_(sixDoFSolver::New(dict.subDict("solver"), *this))
{
addRestraints(dict);
@ -178,6 +180,7 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
aRelax_(sDoFRBM.aRelax_),
aDamp_(sDoFRBM.aDamp_),
report_(sDoFRBM.report_),
updateConstraints_(sDoFRBM.updateConstraints_),
solver_(sDoFRBM.solver_.clone())
{}
@ -304,6 +307,31 @@ void Foam::sixDoFRigidBodyMotion::updateAcceleration
}
void Foam::sixDoFRigidBodyMotion::updateConstraints()
{
if (!updateConstraints_)
{
return;
}
pointConstraint pct;
pointConstraint pcr;
forAll(constraints_, i)
{
constraints_[i].setCentreOfRotation(initialCentreOfRotation_);
constraints_[i].constrainTranslation(pct);
constraints_[i].constrainRotation(pcr);
}
tConstraints_ = pct.constraintTransformation();
rConstraints_ = pcr.constraintTransformation();
Info<< "Translational constraint tensor " << tConstraints_ << nl
<< "Rotational constraint tensor " << rConstraints_ << endl;
}
void Foam::sixDoFRigidBodyMotion::update
(
bool firstIter,

11
src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2015 OpenFOAM Foundation
Copyright (C) 2016-2019 OpenCFD Ltd.
Copyright (C) 2016-2023 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -119,6 +119,9 @@ class sixDoFRigidBodyMotion
//- Reporting of motion data (on or off)
bool report_;
//- Flag to enable time-variant constraints
bool updateConstraints_;
//- Motion solver
autoPtr<sixDoFSolver> solver_;
@ -152,6 +155,9 @@ class sixDoFRigidBodyMotion
//- Update and relax accelerations from the force and torque
void updateAcceleration(const vector& fGlobal, const vector& tauGlobal);
//- Update the constraints to the object
void updateConstraints();
// Access functions retained as private because of the risk of
// confusion over what is a body local frame vector and what is global
@ -271,6 +277,9 @@ public:
//- Return the report Switch
inline bool report() const;
//- Return the update-constraints flag
inline bool updateConstraints() const;
//- Return time
inline const Time& time() const;

7
src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2015 OpenFOAM Foundation
Copyright (C) 2016-2019 OpenCFD Ltd.
Copyright (C) 2016-2023 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -275,6 +275,11 @@ inline bool Foam::sixDoFRigidBodyMotion::report() const
return report_;
}
inline bool Foam::sixDoFRigidBodyMotion::updateConstraints() const
{
return updateConstraints_;
}
inline void Foam::sixDoFRigidBodyMotion::newTime()
{

4
src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionIO.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2014 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
Copyright (C) 2020-2023 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -39,6 +39,7 @@ bool Foam::sixDoFRigidBodyMotion::read(const dictionary& dict)
aRelax_ = dict.getOrDefault<scalar>("accelerationRelaxation", 1);
aDamp_ = dict.getOrDefault<scalar>("accelerationDamping", 1);
report_ = dict.getOrDefault<Switch>("report", false);
updateConstraints_ = dict.getOrDefault("updateConstraints", false);
restraints_.clear();
addRestraints(dict);
@ -61,6 +62,7 @@ void Foam::sixDoFRigidBodyMotion::write(Ostream& os) const
os.writeEntry("accelerationRelaxation", aRelax_);
os.writeEntry("accelerationDamping", aDamp_);
os.writeEntry("report", report_);
os.writeEntry("updateConstraints", updateConstraints_);
if (!restraints_.empty())
{

5
src/sixDoFRigidBodyMotion/sixDoFSolvers/CrankNicolson/CrankNicolson.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2015 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
Copyright (C) 2020-2023 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -75,6 +75,9 @@ void Foam::sixDoFSolvers::CrankNicolson::solve
// Update the linear acceleration and torque
updateAcceleration(fGlobal, tauGlobal);
// Update the constraints to the object
updateConstraints();
// Correct linear velocity
v() = tConstraints()
& (v0() + aDamp()*deltaT*(aoc_*a() + (1 - aoc_)*a0()));

5
src/sixDoFRigidBodyMotion/sixDoFSolvers/Newmark/Newmark.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2015-2016 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
Copyright (C) 2020-2023 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -82,6 +82,9 @@ void Foam::sixDoFSolvers::Newmark::solve
// Update the linear acceleration and torque
updateAcceleration(fGlobal, tauGlobal);
// Update the constraints to the object
updateConstraints();
// Correct linear velocity
v() =
tConstraints()

4
src/sixDoFRigidBodyMotion/sixDoFSolvers/sixDoFSolver/sixDoFSolver.H
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2015 OpenFOAM Foundation
Copyright (C) 2023 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -128,6 +129,9 @@ protected:
const vector& tauGlobal
);
//- Update the constraints to the object
inline void updateConstraints();
public:

6
src/sixDoFRigidBodyMotion/sixDoFSolvers/sixDoFSolver/sixDoFSolverI.H
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2015 OpenFOAM Foundation
Copyright (C) 2023 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -129,5 +130,10 @@ inline void Foam::sixDoFSolver::updateAcceleration
body_.updateAcceleration(fGlobal, tauGlobal);
}
inline void Foam::sixDoFSolver::updateConstraints()
{
body_.updateConstraints();
}
// ************************************************************************* //

4
src/sixDoFRigidBodyMotion/sixDoFSolvers/symplectic/symplectic.C
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2015-2016 OpenFOAM Foundation
Copyright (C) 2023 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -85,6 +86,9 @@ void Foam::sixDoFSolvers::symplectic::solve
// Update the linear acceleration and torque
updateAcceleration(fGlobal, tauGlobal);
// Update the constraints to the object
updateConstraints();
// Second simplectic step:
// Complete update of linear and angular velocities