Commit f4202d9e authored by Henry Weller's avatar Henry Weller
Browse files

sixDoFSolver: Run-time selectable solver (integrator) for sixDoFRigidBodyMotion

The built-in explicit symplectic integrator has been replaced by a
general framework supporting run-time selectable integrators.  Currently
the explicit symplectic, implicit Crank-Nicolson and implicit Newmark
methods are provided, all of which are 2nd-order in time:

Symplectic 2nd-order explicit time-integrator for 6DoF solid-body motion:

    Reference:
        Dullweber, A., Leimkuhler, B., & McLachlan, R. (1997).
        Symplectic splitting methods for rigid body molecular dynamics.
        The Journal of chemical physics, 107(15), 5840-5851.

    Can only be used for explicit integration of the motion of the body,
    i.e. may only be called once per time-step, no outer-correctors may be
    applied.  For implicit integration with outer-correctors choose either
    CrankNicolson or Newmark schemes.

    Example specification in dynamicMeshDict:
    solver
    {
        type    symplectic;
    }

Newmark 2nd-order time-integrator for 6DoF solid-body motion:

    Reference:
        Newmark, N. M. (1959).
        A method of computation for structural dynamics.
        Journal of the Engineering Mechanics Division, 85(3), 67-94.

    Example specification in dynamicMeshDict:
    solver
    {
        type    Newmark;
        gamma   0.5;    // Velocity integration coefficient
        beta    0.25;   // Position integration coefficient
    }

Crank-Nicolson 2nd-order time-integrator for 6DoF solid-body motion:

    The off-centering coefficients for acceleration (velocity integration) and
    velocity (position/orientation integration) may be specified but default
    values of 0.5 for each are used if they are not specified.  With the default
    off-centering this scheme is equivalent to the Newmark scheme with default
    coefficients.

    Example specification in dynamicMeshDict:
    solver
    {
        type    CrankNicolson;
        aoc     0.5;    // Acceleration off-centering coefficient
        voc     0.5;    // Velocity off-centering coefficient
    }

Both the Newmark and Crank-Nicolson are proving more robust and reliable
than the symplectic method for solving complex coupled problems and the
tutorial cases have been updated to utilize this.

In this new framework it would be straight forward to add other methods
should the need arise.

Henry G. Weller
CFD Direct
parent 11f35cc7
......@@ -133,7 +133,11 @@ int main(int argc, char *argv[])
{
Info<< "Time = " << runTime.timeName() << endl;
mesh.update();
for (int i = 0; i<2; i++)
{
mesh.update();
}
mesh.checkMesh(true);
if (checkAMI)
......
......@@ -31,4 +31,10 @@ sixDoFRigidBodyMotionSolver/sixDoFRigidBodyMotionSolver.C
sixDoFRigidBodyMotionSolver/externalPointEdgePoint.C
sixDoFRigidBodyMotionSolver/pointPatchDist.C
sixDoFSolvers/sixDoFSolver/sixDoFSolver.C
sixDoFSolvers/sixDoFSolver/newSixDoFSolver.C
sixDoFSolvers/symplectic/symplectic.C
sixDoFSolvers/CrankNicolson/CrankNicolson.C
sixDoFSolvers/Newmark/Newmark.C
LIB = $(FOAM_LIBBIN)/libsixDoFRigidBodyMotion
......@@ -224,11 +224,6 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
f.calcForcesMoment();
// Patch force data is valid for the current positions, so
// calculate the forces on the motion object from this data, then
// update the positions
motion_.updatePosition(firstIter, t.deltaTValue(), t.deltaT0Value());
// Get the forces on the patch faces at the current positions
if (lookupGravity_ == 1)
......@@ -242,7 +237,7 @@ void sixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
// scalar ramp = min(max((t.value() - 5)/10, 0), 1);
scalar ramp = 1.0;
motion_.updateAcceleration
motion_.update
(
firstIter,
ramp*(f.forceEff() + motion_.mass()*g_),
......
......@@ -160,8 +160,6 @@ void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
firstIter = true;
}
motion_.updatePosition(firstIter, t.deltaTValue(), t.deltaT0Value());
vector gravity = vector::zero;
if (db().foundObject<uniformDimensionedVectorField>("g"))
......@@ -173,7 +171,7 @@ void uncoupledSixDoFRigidBodyDisplacementPointPatchVectorField::updateCoeffs()
}
// Do not modify the accelerations, except with gravity, where available
motion_.updateAcceleration
motion_.update
(
firstIter,
gravity*motion_.mass(),
......
......@@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "sixDoFRigidBodyMotion.H"
#include "sixDoFSolver.H"
#include "septernion.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
......@@ -84,7 +85,8 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion()
momentOfInertia_(diagTensor::one*VSMALL),
aRelax_(1.0),
aDamp_(1.0),
report_(false)
report_(false),
solver_(NULL)
{}
......@@ -121,7 +123,8 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
momentOfInertia_(dict.lookup("momentOfInertia")),
aRelax_(dict.lookupOrDefault<scalar>("accelerationRelaxation", 1.0)),
aDamp_(dict.lookupOrDefault<scalar>("accelerationDamping", 1.0)),
report_(dict.lookupOrDefault<Switch>("report", false))
report_(dict.lookupOrDefault<Switch>("report", false)),
solver_(sixDoFSolver::New(dict.subDict("solver"), *this))
{
addRestraints(dict);
......@@ -171,6 +174,12 @@ Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::sixDoFRigidBodyMotion::~sixDoFRigidBodyMotion()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::sixDoFRigidBodyMotion::addRestraints
......@@ -257,50 +266,35 @@ void Foam::sixDoFRigidBodyMotion::addConstraints
}
void Foam::sixDoFRigidBodyMotion::updatePosition
void Foam::sixDoFRigidBodyMotion::updateAcceleration
(
bool firstIter,
scalar deltaT,
scalar deltaT0
const vector& fGlobal,
const vector& tauGlobal
)
{
if (Pstream::master())
{
if (firstIter)
{
// First simplectic step:
// Half-step for linear and angular velocities
// Update position and orientation
v() = tConstraints_ & (v0() + aDamp_*0.5*deltaT0*a());
pi() = rConstraints_ & (pi0() + aDamp_*0.5*deltaT0*tau());
centreOfRotation() = centreOfRotation0() + deltaT*v();
}
else
{
// For subsequent iterations use Crank-Nicolson
static bool first = false;
v() = tConstraints_
& (v0() + aDamp_*0.5*deltaT*(a() + motionState0_.a()));
pi() = rConstraints_
& (pi0() + aDamp_*0.5*deltaT*(tau() + motionState0_.tau()));
// Save the previous iteration accelerations for relaxation
vector aPrevIter = a();
vector tauPrevIter = tau();
centreOfRotation() =
centreOfRotation0() + 0.5*deltaT*(v() + motionState0_.v());
}
// Calculate new accelerations
a() = fGlobal/mass_;
tau() = (Q().T() & tauGlobal);
applyRestraints();
// Correct orientation
Tuple2<tensor, vector> Qpi = rotate(Q0(), pi(), deltaT);
Q() = Qpi.first();
pi() = rConstraints_ & Qpi.second();
// Relax accelerations on all but first iteration
if (!first)
{
a() = aRelax_*a() + (1 - aRelax_)*aPrevIter;
tau() = aRelax_*tau() + (1 - aRelax_)*tauPrevIter;
}
Pstream::scatter(motionState_);
first = false;
}
void Foam::sixDoFRigidBodyMotion::updateAcceleration
void Foam::sixDoFRigidBodyMotion::update
(
bool firstIter,
const vector& fGlobal,
......@@ -309,44 +303,9 @@ void Foam::sixDoFRigidBodyMotion::updateAcceleration
scalar deltaT0
)
{
static bool first = false;
if (Pstream::master())
{
// Save the previous iteration accelerations for relaxation
vector aPrevIter = a();
vector tauPrevIter = tau();
// Calculate new accelerations
a() = fGlobal/mass_;
tau() = (Q().T() & tauGlobal);
applyRestraints();
// Relax accelerations on all but first iteration
if (!first)
{
a() = aRelax_*a() + (1 - aRelax_)*aPrevIter;
tau() = aRelax_*tau() + (1 - aRelax_)*tauPrevIter;
}
first = false;
if (firstIter)
{
// Second simplectic step:
// Complete update of linear and angular velocities
v() += tConstraints_ & aDamp_*0.5*deltaT*a();
pi() += rConstraints_ & aDamp_*0.5*deltaT*tau();
}
else
{
// For subsequent iterations use Crank-Nicolson
v() = tConstraints_
& (v0() + aDamp_*0.5*deltaT*(a() + motionState0_.a()));
pi() = rConstraints_
& (pi0() + aDamp_*0.5*deltaT*(tau() + motionState0_.tau()));
}
solver_->solve(firstIter, fGlobal, tauGlobal, deltaT, deltaT0);
if (report_)
{
......@@ -358,7 +317,6 @@ void Foam::sixDoFRigidBodyMotion::updateAcceleration
}
void Foam::sixDoFRigidBodyMotion::status() const
{
Info<< "6-DoF rigid body motion" << nl
......
......@@ -25,28 +25,16 @@ Class
Foam::sixDoFRigidBodyMotion
Description
Six degree of freedom motion for a rigid body. Angular momentum
stored in body fixed reference frame. Reference orientation of
the body (where Q = I) must align with the cartesian axes such
that the Inertia tensor is in principle component form.
Six degree of freedom motion for a rigid body.
Symplectic motion as per:
Angular momentum stored in body fixed reference frame. Reference
orientation of the body (where Q = I) must align with the cartesian axes
such that the Inertia tensor is in principle component form. Can add
restraints (e.g. a spring) and constraints (e.g. motion may only be on a
plane).
title = {Symplectic splitting methods for rigid body molecular dynamics},
publisher = {AIP},
year = {1997},
journal = {The Journal of Chemical Physics},
volume = {107},
number = {15},
pages = {5840-5851},
url = {http://link.aip.org/link/?JCP/107/5840/1},
doi = {10.1063/1.474310}
Changes to Crank-Nicolson integration if subsequent iterations to handle
implicit coupling between pressure and motion.
Can add restraints (e.g. a spring)
and constraints (e.g. motion may only be on a plane).
The time-integrator for the motion is run-time selectable with options for
symplectic (explicit), Crank-Nicolson and Newmark schemes.
SourceFiles
sixDoFRigidBodyMotionI.H
......@@ -69,12 +57,17 @@ SourceFiles
namespace Foam
{
// Forward declarations
class sixDoFSolver;
/*---------------------------------------------------------------------------*\
Class sixDoFRigidBodyMotion Declaration
\*---------------------------------------------------------------------------*/
class sixDoFRigidBodyMotion
{
friend class sixDoFSolver;
// Private data
//- Motion state data object
......@@ -120,6 +113,9 @@ class sixDoFRigidBodyMotion
//- Switch to turn reporting of motion data on and off
Switch report_;
//- Motion solver
autoPtr<sixDoFSolver> solver_;
// Private Member Functions
......@@ -147,6 +143,9 @@ class sixDoFRigidBodyMotion
//- Apply the restraints to the object
void applyRestraints();
//- Update and relax accelerations from the force and torque
void updateAcceleration(const vector& fGlobal, const vector& tauGlobal);
// Access functions retained as private because of the risk of
// confusion over what is a body local frame vector and what is global
......@@ -179,21 +178,6 @@ class sixDoFRigidBodyMotion
//- Return the current torque
inline const vector& tau() const;
//- Return the orientation at previous time-step
inline const tensor& Q0() const;
//- Return the velocity at previous time-step
inline const vector& v0() const;
//- Return the acceleration at previous time-step
inline const vector& a0() const;
//- Return the angular momentum at previous time-step
inline const vector& pi0() const;
//- Return the torque at previous time-step
inline const vector& tau0() const;
// Edit
......@@ -237,6 +221,10 @@ public:
sixDoFRigidBodyMotion(const sixDoFRigidBodyMotion&);
//- Destructor
~sixDoFRigidBodyMotion();
// Member Functions
// Access
......@@ -253,9 +241,6 @@ public:
//- Return the current centre of rotation
inline const point& centreOfRotation() const;
//- Return the centre of rotation at previous time-step
inline const point& centreOfRotation0() const;
//- Return the initial centre of mass
inline const point& initialCentreOfMass() const;
......@@ -301,16 +286,9 @@ public:
// Update state
//- First symplectic velocity adjust and motion part, required
// before force calculation. Takes old timestep for variable
// timestep cases.
// Changes to Crank-Nicolson integration for subsequent iterations.
void updatePosition(bool firstIter, scalar deltaT, scalar deltaT0);
//- Second symplectic velocity adjust part
// required after motion and force calculation.
//- Symplectic integration of velocities, orientation and position.
// Changes to Crank-Nicolson integration for subsequent iterations.
void updateAcceleration
void update
(
bool firstIter,
const vector& fGlobal,
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
......@@ -161,36 +161,6 @@ inline const Foam::vector& Foam::sixDoFRigidBodyMotion::tau() const
}
inline const Foam::tensor& Foam::sixDoFRigidBodyMotion::Q0() const
{
return motionState0_.Q();
}
inline const Foam::vector& Foam::sixDoFRigidBodyMotion::v0() const
{
return motionState0_.v();
}
inline const Foam::vector& Foam::sixDoFRigidBodyMotion::a0() const
{
return motionState0_.a();
}
inline const Foam::vector& Foam::sixDoFRigidBodyMotion::pi0() const
{
return motionState0_.pi();
}
inline const Foam::vector& Foam::sixDoFRigidBodyMotion::tau0() const
{
return motionState0_.tau();
}
inline Foam::point& Foam::sixDoFRigidBodyMotion::initialCentreOfRotation()
{
return initialCentreOfRotation_;
......@@ -261,12 +231,6 @@ inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfRotation() const
}
inline const Foam::point& Foam::sixDoFRigidBodyMotion::centreOfRotation0() const
{
return motionState0_.centreOfRotation();
}
inline const Foam::point&
Foam::sixDoFRigidBodyMotion::initialCentreOfMass() const
{
......
......@@ -205,9 +205,7 @@ void Foam::sixDoFRigidBodyMotionSolver::solve()
if (test_)
{
motion_.updatePosition(firstIter, t.deltaTValue(), t.deltaT0Value());
motion_.updateAcceleration
motion_.update
(
firstIter,
ramp*(motion_.mass()*g.value()),
......@@ -230,12 +228,7 @@ void Foam::sixDoFRigidBodyMotionSolver::solve()
f.calcForcesMoment();
// Patch force data is valid for the current positions, so
// calculate the forces on the motion object from this data, then
// update the positions
motion_.updatePosition(firstIter, t.deltaTValue(), t.deltaT0Value());
motion_.updateAcceleration
motion_.update
(
firstIter,
ramp*(f.forceEff() + motion_.mass()*g.value()),
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM 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.
OpenFOAM 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 OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "CrankNicolson.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
namespace sixDoFSolvers
{
defineTypeNameAndDebug(CrankNicolson, 0);
addToRunTimeSelectionTable(sixDoFSolver, CrankNicolson, dictionary);
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::sixDoFSolvers::CrankNicolson::CrankNicolson
(
const dictionary& dict,
sixDoFRigidBodyMotion& body
)
:
sixDoFSolver(body),
aoc_(dict.lookupOrDefault<scalar>("aoc", 0.5)),
voc_(dict.lookupOrDefault<scalar>("voc", 0.5))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::sixDoFSolvers::CrankNicolson::~CrankNicolson()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::sixDoFSolvers::CrankNicolson::solve
(
bool firstIter,
const vector& fGlobal,
const vector& tauGlobal,
scalar deltaT,
scalar deltaT0
)
{
// Update the linear acceleration and torque
updateAcceleration(fGlobal, tauGlobal);
// Correct linear velocity
v() = tConstraints()
& (v0() + aDamp()*deltaT*(aoc_*a() + (1 - aoc_)*a0()));
// Correct angular momentum
pi() = rConstraints()
& (pi0() + aDamp()*deltaT*(aoc_*tau() + (1 - aoc_)*tau0()));
// Correct position
centreOfRotation() =
centreOfRotation0() + deltaT*(voc_*v() + (1 - voc_)*v0());
// Correct orientation
Tuple2<tensor, vector> Qpi =
rotate(Q0(), (voc_*pi() + (1 - voc_)*pi0()), deltaT);
Q() = Qpi.first();
}
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM 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.
OpenFOAM 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 OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::sixDoFSolvers::CrankNicolson
Description
Crank-Nicolson 2nd-order time-integrator for 6DoF solid-body motion.
The off-centering coefficients for acceleration (velocity integration) and
velocity (position/orientation integration) may be specified but default
values of 0.5 for each are used if they are not specified. With the default
off-centering this scheme is equivalent to the Newmark scheme with default
coefficients.
Example specification in dynamicMeshDict:
\verbatim
solver
{
type CrankNicolson;
aoc 0.5; // Acceleration off-centering coefficient
voc 0.5; // Velocity off-centering coefficient
}
\endverbatim
SeeAlso
Foam::sixDoFSolvers::Newmark
SourceFiles
CrankNicolson.C
\*---------------------------------------------------------------------------*/