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

rigidBodyDynamics/rigidBodySolvers: Added run-time selectable solvers to...

rigidBodyDynamics/rigidBodySolvers: Added run-time selectable solvers to integrate the rigid-body motion

Currently supported solvers: symplectic, Newmark, CrankNicolson

The symplectic solver should only be used if iteration over the forces
and body-motion is not required.  Newmark and CrankNicolson both require
iteration to provide 2nd-order behavior.

See applications/test/rigidBodyDynamics/spring for an example of the
application of the Newmark solver.

This development is sponsored by Carnegie Wave Energy Ltd.
parent f4baa3a8
solver
{
type Newmark;
}
// It is necessary to iterate for the Newmark solver
nIter 2;
bodies
{
weight
......
......@@ -29,7 +29,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "rigidBodyModel.H"
#include "rigidBodyMotion.H"
#include "masslessBody.H"
#include "sphere.H"
#include "joints.H"
......@@ -45,18 +45,19 @@ using namespace RBD;
int main(int argc, char *argv[])
{
dictionary springDict(IFstream("spring")());
// Create the spring model from dictionary
rigidBodyModel spring(dictionary(IFstream("spring")()));
rigidBodyMotion spring(springDict);
Info<< spring << endl;
label nIter(readLabel(springDict.lookup("nIter")));
// Create the joint-space state fields
rigidBodyModelState springState(spring);
scalarField& q = springState.q();
scalarField& qDot = springState.qDot();
scalarField& qDdot = springState.qDdot();
Info<< spring << endl;
// Create the joint-space force field
scalarField tau(spring.nDoF(), Zero);
// Create the external body force field
Field<spatialVector> fx(spring.nBodies(), Zero);
OFstream qFile("qVsTime");
......@@ -66,27 +67,23 @@ int main(int argc, char *argv[])
scalar deltaT = 0.002;
for (scalar t=0; t<4; t+=deltaT)
{
qDot += 0.5*deltaT*qDdot;
q += deltaT*qDot;
// Update the body-state prior to the evaluation of the restraints
spring.forwardDynamicsCorrection(springState);
// Accumulate the restraint forces
fx = Zero;
spring.applyRestraints(fx);
// Calculate the body acceleration for the given state
// and restraint forces
spring.forwardDynamics(springState, tau, fx);
// Update the velocity
qDot += 0.5*deltaT*qDdot;
spring.newTime();
for (label i=0; i<nIter; i++)
{
spring.update
(
deltaT,
deltaT,
tau,
fx
);
}
// Write the results for graph generation
// using 'gnuplot spring.gnuplot'
qFile << t << " " << q[0] << endl;
qDotFile << t << " " << qDot[0] << endl;
qFile << t << " " << spring.state().q()[0] << endl;
qDotFile << t << " " << spring.state().qDot()[0] << endl;
}
Info<< "\nEnd\n" << endl;
......
......@@ -38,4 +38,13 @@ rigidBodyModel/forwardDynamics.C
rigidBodyModelState/rigidBodyModelState.C
rigidBodyModelState/rigidBodyModelStateIO.C
rigidBodyMotion/rigidBodyMotion.C
rigidBodyMotion/rigidBodyMotionIO.C
rigidBodySolvers/rigidBodySolver/rigidBodySolver.C
rigidBodySolvers/rigidBodySolver/newRigidBodySolver.C
rigidBodySolvers/symplectic/symplectic.C
rigidBodySolvers/Newmark/Newmark.C
rigidBodySolvers/CrankNicolson/CrankNicolson.C
LIB = $(FOAM_LIBBIN)/librigidBodyDynamics
......@@ -22,7 +22,7 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::rigidBodyModelState
Foam::RBD::rigidBodyModelState
Description
Holds the motion state of rigid-body model.
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 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 "rigidBodyMotion.H"
#include "rigidBodySolver.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::RBD::rigidBodyMotion::rigidBodyMotion()
:
rigidBodyModel(),
motionState_(*this),
motionState0_(*this),
aRelax_(1.0),
aDamp_(1.0),
report_(false),
solver_(NULL)
{}
Foam::RBD::rigidBodyMotion::rigidBodyMotion
(
const dictionary& dict
)
:
rigidBodyModel(dict),
motionState_(*this),
motionState0_(motionState_),
aRelax_(dict.lookupOrDefault<scalar>("accelerationRelaxation", 1.0)),
aDamp_(dict.lookupOrDefault<scalar>("accelerationDamping", 1.0)),
report_(dict.lookupOrDefault<Switch>("report", false)),
solver_(rigidBodySolver::New(*this, dict.subDict("solver")))
{}
Foam::RBD::rigidBodyMotion::rigidBodyMotion
(
const dictionary& dict,
const dictionary& stateDict
)
:
rigidBodyModel(dict),
motionState_(*this, stateDict),
motionState0_(motionState_),
aRelax_(dict.lookupOrDefault<scalar>("accelerationRelaxation", 1.0)),
aDamp_(dict.lookupOrDefault<scalar>("accelerationDamping", 1.0)),
report_(dict.lookupOrDefault<Switch>("report", false)),
solver_(rigidBodySolver::New(*this, dict.subDict("solver")))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::RBD::rigidBodyMotion::~rigidBodyMotion()
{}
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
void Foam::RBD::rigidBodyMotion::update
(
scalar deltaT,
scalar deltaT0,
const scalarField& tau,
const Field<spatialVector>& fx
)
{
if (Pstream::master())
{
solver_->solve(deltaT, deltaT0, tau, fx);
if (report_)
{
status();
}
}
Pstream::scatter(motionState_);
}
void Foam::RBD::rigidBodyMotion::status() const
{
/*
Info<< "Rigid body motion" << nl
<< " Centre of rotation: " << centreOfRotation() << nl
<< " Centre of mass: " << centreOfMass() << nl
<< " Orientation: " << orientation() << nl
<< " Linear velocity: " << v() << nl
<< " Angular velocity: " << omega()
<< endl;
*/
}
/*
Foam::tmp<Foam::pointField> Foam::RBD::rigidBodyMotion::transform
(
const pointField& initialPoints
) const
{
return
(
centreOfRotation()
+ (Q() & initialQ_.T() & (initialPoints - initialCentreOfRotation_))
);
}
Foam::tmp<Foam::pointField> Foam::RBD::rigidBodyMotion::transform
(
const pointField& initialPoints,
const scalarField& scale
) const
{
// Calculate the transformation septerion from the initial state
septernion s
(
centreOfRotation() - initialCentreOfRotation(),
quaternion(Q() & initialQ().T())
);
tmp<pointField> tpoints(new pointField(initialPoints));
pointField& points = tpoints.ref();
forAll(points, pointi)
{
// Move non-stationary points
if (scale[pointi] > SMALL)
{
// Use solid-body motion where scale = 1
if (scale[pointi] > 1 - SMALL)
{
points[pointi] = transform(initialPoints[pointi]);
}
// Slerp septernion interpolation
else
{
septernion ss(slerp(septernion::I, s, scale[pointi]));
points[pointi] =
initialCentreOfRotation()
+ ss.transform
(
initialPoints[pointi]
- initialCentreOfRotation()
);
}
}
}
return tpoints;
}
*/
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 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::RBD::rigidBodyMotion
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. 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
rigidBodyMotionI.H
rigidBodyMotion.C
rigidBodyMotionIO.C
\*---------------------------------------------------------------------------*/
#ifndef rigidBodyMotion_H
#define rigidBodyMotion_H
#include "rigidBodyModel.H"
#include "rigidBodyModelState.H"
#include "pointField.H"
#include "Switch.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace RBD
{
// Forward declarations
class rigidBodySolver;
/*---------------------------------------------------------------------------*\
Class rigidBodyMotion Declaration
\*---------------------------------------------------------------------------*/
class rigidBodyMotion
:
public rigidBodyModel
{
friend class rigidBodySolver;
// Private data
//- Motion state data object
rigidBodyModelState motionState_;
//- Motion state data object for previous time-step
rigidBodyModelState motionState0_;
//- Acceleration relaxation coefficient
scalar aRelax_;
//- Acceleration damping coefficient (for steady-state simulations)
scalar aDamp_;
//- Switch to turn reporting of motion data on and off
Switch report_;
//- Motion solver
autoPtr<rigidBodySolver> solver_;
//- Construct as copy
rigidBodyMotion(const rigidBodyMotion&);
public:
// Constructors
//- Construct null
rigidBodyMotion();
//- Construct from dictionary
rigidBodyMotion
(
const dictionary& dict
);
//- Construct from constant and state dictionaries
rigidBodyMotion
(
const dictionary& dict,
const dictionary& stateDict
);
//- Destructor
~rigidBodyMotion();
// Member Functions
// Access
//- Return the report Switch
inline bool report() const;
//- Return the motion state
inline const rigidBodyModelState& state() const;
// Edit
//- Store the motion state at the beginning of the time-step
inline void newTime();
// Update state
//- Integration of velocities, orientation and position.
void update
(
scalar deltaT,
scalar deltaT0,
const scalarField& tau,
const Field<spatialVector>& fx
);
//- Report the status of the motion
void status() const;
// Transformations
/*
//- Transform the given initial state point by the current motion
// state
inline point transform(const point& initialPoints) const;
//- Transform the given initial state pointField by the current
// motion state
tmp<pointField> transform(const pointField& initialPoints) const;
//- Transform the given initial state pointField by the current
// motion state scaled by the given scale
tmp<pointField> transform
(
const pointField& initialPoints,
const scalarField& scale
) const;
*/
//- Write
void write(Ostream&) const;
//- Read coefficients dictionary and update system parameters,
// constraints and restraints but not the current state
bool read(const dictionary& dict);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace RBD
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "rigidBodyMotionI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 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/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
inline bool Foam::RBD::rigidBodyMotion::report() const
{
return report_;
}
inline const Foam::RBD::rigidBodyModelState&
Foam::RBD::rigidBodyMotion::state() const
{
return motionState_;
}
inline void Foam::RBD::rigidBodyMotion::newTime()
{
motionState0_ = motionState_;
}
/*
inline Foam::point Foam::RBD::rigidBodyMotion::transform
(
const point& initialPoint
) const
{
return
(
centreOfRotation()
+ (Q() & initialQ_.T() & (initialPoint - initialCentreOfRotation_))
);
}
*/
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2016 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 "rigidBodyMotion.H"
#include "IOstreams.H"
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
bool Foam::RBD::rigidBodyMotion::read(const dictionary& dict)
{
rigidBodyModel::read(dict);
aRelax_ = dict.lookupOrDefault<scalar>("accelerationRelaxation", 1.0);
aDamp_ = dict.lookupOrDefault<scalar>("accelerationDamping", 1.0);
report_ = dict.lookupOrDefault<Switch>("report", false);
return true;
}
void Foam::RBD::rigidBodyMotion::write(Ostream& os) const
{
rigidBodyModel::write(os);
os.writeKeyword("accelerationRelaxation")
<< aRelax_ << token::END_STATEMENT << nl;
os.writeKeyword("accelerationDamping")
<< aDamp_ << token::END_STATEMENT << nl;
os.writeKeyword("report")
<< report_ << token::END_STATEMENT << nl;
}
// ************************************************************************* //