diff --git a/src/ODE/Make/files b/src/ODE/Make/files index 9c2455c38..dc7197423 100644 --- a/src/ODE/Make/files +++ b/src/ODE/Make/files @@ -9,12 +9,12 @@ $(ODESolvers)/SIBS/SIBS.C $(ODESolvers)/SIBS/SIMPR.C $(ODESolvers)/SIBS/polyExtrapolate.C -translationODE = translationODE -$(translationODE)/translationODE.C - sixDOF = sixDOF $(sixDOF)/finiteRotation/finiteRotation.C $(sixDOF)/sixDOFqODE/sixDOFqODE.C $(sixDOF)/sixDOFbodies/sixDOFbodies.C +$(sixDOF)/sixDOFODE/sixDOFODE.C +$(sixDOF)/sixDOFODE/newSixDOFODE.C + LIB = $(FOAM_LIBBIN)/libODE diff --git a/src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C b/src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C new file mode 100644 index 000000000..9255dc910 --- /dev/null +++ b/src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C @@ -0,0 +1,73 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright held by original author + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "error.H" +#include "objectRegistry.H" +#include "sixDOFODE.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +Foam::autoPtr Foam::sixDOFODE::New(const IOobject& io) +{ + word sixDOFODETypeName; + + // Get object registry + const objectRegistry& database = io.db(); + + // Check whether the dictionary is in the database + if (database.foundObject(io.name())) + { + sixDOFODETypeName = + word + ( + database.lookupObject(io.name()).lookup("type") + ); + } + else + { + sixDOFODETypeName = word(IOdictionary(io).lookup("type")); + } + + dictionaryConstructorTable::iterator cstrIter = + dictionaryConstructorTablePtr_->find(sixDOFODETypeName); + + if (cstrIter == dictionaryConstructorTablePtr_->end()) + { + FatalErrorIn + ( + "sixDOFODE::New(const IOobject& io)" + ) << "Unknown sixDOFODE " << sixDOFODETypeName + << endl << endl + << "Valid sixDOFODE types are:" << endl + << dictionaryConstructorTablePtr_->toc() + << exit(FatalError); + } + + return autoPtr(cstrIter()(io)); +} + + +// ************************************************************************* // diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C new file mode 100644 index 000000000..1df3ce5df --- /dev/null +++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.C @@ -0,0 +1,252 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Class + sixDOFODE + +Description + Abstract base class for six-degrees-of-freedom (6DOF) ordinary differential + equations + +Author + Dubravko Matijasevic, FSB Zagreb. All rights reserved. + Hrvoje Jasak, FSB Zagreb. All rights reserved. + Vuko Vukcevic, FSB Zagreb. All rights reserved. + +\*---------------------------------------------------------------------------*/ + +#include "objectRegistry.H" +#include "sixDOFODE.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +defineTypeNameAndDebug(sixDOFODE, 0); +defineRunTimeSelectionTable(sixDOFODE, dictionary); + +} + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + +void Foam::sixDOFODE::aitkensRelaxation +( + const scalar min, + const scalar max +) +{ + // Calculate translational relax factor + const scalar saveOldRelFacT = oldRelaxFactorT_; + oldRelaxFactorT_ = relaxFactorT_; + + if(magSqr(A_[0] - An_[1] - A_[1] - An_[2]) > SMALL) + { + relaxFactorT_ = saveOldRelFacT + (saveOldRelFacT - 1)* + ( + (A_[1] - An_[2]) + & (A_[0] - An_[1] - A_[1] - An_[2]) + )/ + magSqr(A_[0] - An_[1] - A_[1] - An_[2]); + } + else + { + relaxFactorT_ = min; + } + + // Bound the relaxation factor for stability + if (relaxFactorT_ > max) + { + relaxFactorT_ = max; + } + else if (relaxFactorT_ < min) + { + relaxFactorT_ = min; + } + + // Calculate rotational relax factor + const scalar saveOldRelFacR = oldRelaxFactorR_; + oldRelaxFactorR_ = relaxFactorR_; + + if + ( + magSqr(OmegaDot_[0] - OmegaDotn_[1] - OmegaDot_[1] - OmegaDotn_[2]) + > SMALL + ) + { + relaxFactorR_ = saveOldRelFacR + + (saveOldRelFacR - 1)* + ( + (OmegaDot_[1] - OmegaDotn_[2]) + & (OmegaDot_[0] - OmegaDotn_[1] - OmegaDot_[1] - OmegaDotn_[2]) + )/ + magSqr(OmegaDot_[0] - OmegaDotn_[1] - OmegaDot_[1] - OmegaDotn_[2]); + } + else + { + relaxFactorR_ = min; + } + + // Bound the relaxation factor for stability + if(relaxFactorR_ > max) + { + relaxFactorR_ = max; + } + else if(relaxFactorR_ < min) + { + relaxFactorR_ = min; + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::sixDOFODE::sixDOFODE(const IOobject& io) +: + IOdictionary(io), + + mass_(lookup("mass")), + momentOfInertia_(lookup("momentOfInertia")), + + Xequilibrium_(lookup("equilibriumPosition")), + linSpringCoeffs_(lookup("linearSpring")), + linDampingCoeffs_(lookup("linearDamping")), + + relaxFactorT_(1.0), + relaxFactorR_(1.0), + oldRelaxFactorT_(1.0), + oldRelaxFactorR_(1.0), + + A_(3, vector::zero), + OmegaDot_(3, vector::zero), + An_(3, vector::zero), + OmegaDotn_(3, vector::zero), + + force_(lookup("force")), + moment_(lookup("moment")), + forceRelative_(lookup("forceRelative")), + momentRelative_(lookup("momentRelative")) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::sixDOFODE::~sixDOFODE() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + + +void Foam::sixDOFODE::relaxAcceleration +( + const scalar minRelFactor, + const scalar maxRelFactor +) +{ + if (mag(minRelFactor - maxRelFactor) < SMALL) + { + // Fixed relaxation + relaxFactorT_ = minRelFactor; + relaxFactorR_ = minRelFactor; + } + else + { + // Use Aitkens relaxation + + // Update Aitkens relaxation factor + aitkensRelaxation(minRelFactor, maxRelFactor); + + // Update non relaxed accelerations + An_[1] = An_[2]; + An_[2] = (force_/mass_).value(); + OmegaDotn_[1] = OmegaDotn_[2]; + OmegaDotn_[2] = (inv(momentOfInertia_) & moment_).value(); + } + + const vector Aold = A_[2]; + const vector OmegaDotold = OmegaDot_[2]; + + force_.value() = + Aold*mass_.value() + + relaxFactorT_*(force_.value() - Aold*mass_.value()); + + moment_.value() = + (momentOfInertia_.value() & OmegaDotold) + + relaxFactorR_* + ( + moment_.value() + - (momentOfInertia_.value() & OmegaDotold) + ); + + // Update relaxed old accelerations + A_[0] = A_[1]; + A_[1] = A_[2]; + A_[2] = (force_/mass_).value(); + OmegaDot_[0] = OmegaDot_[1]; + OmegaDot_[1] = OmegaDot_[2]; + OmegaDot_[2] = (inv(momentOfInertia_) & moment_).value(); +} + + +// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // + +bool Foam::sixDOFODE::writeData(Ostream& os) const +{ + os.writeKeyword("mass") << tab << mass_ << token::END_STATEMENT << nl; + os.writeKeyword("momentOfInertia") << tab << momentOfInertia_ + << token::END_STATEMENT << nl << nl; + + os.writeKeyword("equilibriumPosition") << tab << Xequilibrium_ + << token::END_STATEMENT << nl; + os.writeKeyword("linearSpring") << tab << linSpringCoeffs_ + << token::END_STATEMENT << nl; + os.writeKeyword("linearDamping") << tab << linDampingCoeffs_ + << token::END_STATEMENT << nl << nl; + + os.writeKeyword("force") << tab << force_ + << token::END_STATEMENT << nl; + os.writeKeyword("moment") << tab << moment_ + << token::END_STATEMENT << nl; + os.writeKeyword("forceRelative") << tab << forceRelative_ + << token::END_STATEMENT << nl; + os.writeKeyword("momentRelative") << tab << momentRelative_ + << token::END_STATEMENT << endl; + + return os.good(); +} + + +// * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * * // + +Foam::Ostream& Foam::operator<<(Ostream& os, const sixDOFODE& sds) +{ + sds.writeData(os); + + os.check("Ostream& operator<<(Ostream&, const sixDOFODE"); + + return os; +} + + +// ************************************************************************* // diff --git a/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H new file mode 100644 index 000000000..a613b8254 --- /dev/null +++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODE.H @@ -0,0 +1,395 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.0 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Class + sixDOFODE + +Description + Abstract base class for six-degrees-of-freedom (6DOF) ordinary differential + equations with a fat interface in order to provide representation of + rotations in quaternions or rotation tensors. + +Author + Dubravko Matijasevic, FSB Zagreb. All rights reserved. + Hrvoje Jasak, FSB Zagreb. All rights reserved. + Vuko Vukcevic, FSB Zagreb. All rights reserved. + +SourceFiles + sixDOFODEI.H + sixDOFODE.C + newSixDOFODE.C + +\*---------------------------------------------------------------------------*/ + +#ifndef sixDOFODE_H +#define sixDOFODE_H + +#include "ODE.H" +#include "IOdictionary.H" +#include "dimensionedTypes.H" +#include "autoPtr.H" +#include "runTimeSelectionTables.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +class finiteRotation; +class HamiltonRodriguezRot; + +/*---------------------------------------------------------------------------*\ + Class sixDOFODE Declaration +\*---------------------------------------------------------------------------*/ + +class sixDOFODE +: + public IOdictionary, + public ODE +{ + // Private data + + // Body data + + //- Mass + dimensionedScalar mass_; + + //- Rotational moment of inertia around centre of mass + // in body (relative) coordinates - aligned with main axes + dimensionedDiagTensor momentOfInertia_; + + + // Equilibrium position and spring/damping coefficients + + //- Spring equilibrium position for translation + dimensionedVector Xequilibrium_; + + //- Linear spring coeffs + dimensionedDiagTensor linSpringCoeffs_; + + //- Linear damping coeffs + dimensionedDiagTensor linDampingCoeffs_; + + + // Aitkens relaxation data members + + //- Translational relaxation factor + scalar relaxFactorT_; + + //- Rotational relaxation factor + scalar relaxFactorR_; + + //- Old translational relaxation factor + scalar oldRelaxFactorT_; + + //- Old rotational relaxation factor + scalar oldRelaxFactorR_; + + //- Accelerations from previous iterations + // A_[2] is the old value, A_[1] old old, A_[0] old old old + List A_; + List OmegaDot_; + + //- Previos iteration non relaxed accelerations + List An_; + List OmegaDotn_; + + + // External forces + + //- Force driving the motion in inertial coord. sys. + dimensionedVector force_; + + //- Moment driving the motion in inertial coord. sys. + dimensionedVector moment_; + + //- Force driving the motion in relative coord. sys. + dimensionedVector forceRelative_; + + //- Moment driving the motion in relative coord. sys. + dimensionedVector momentRelative_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + sixDOFODE(const sixDOFODE&); + + //- Disallow default bitwise assignment + void operator=(const sixDOFODE&); + + +protected: + + // Protected Member Functions + + //- Set ODE coefficients from position and rotation + virtual void setCoeffs() = 0; + + //- Update Aitkens relaxation parameters + void aitkensRelaxation(const scalar min, const scalar max); + + +public: + + //- Run-time type information + TypeName("sixDOFODE"); + + + // Declare run-time constructor selection table + + declareRunTimeSelectionTable + ( + autoPtr, + sixDOFODE, + dictionary, + (const IOobject& io), + (io) + ); + + + // Constructors + + //- Construct from dictionary + sixDOFODE(const IOobject& io); + + //- Return a clone, changing the name + virtual autoPtr clone(const word& name) const = 0; + + + // Selectors + + //- Return autoPtr to the selected sixDOFODE + static autoPtr New(const IOobject& io); + + + // Destructor + + virtual ~sixDOFODE(); + + + // Member Functions + + // Access to common data + + //- Return mass + inline const dimensionedScalar& mass() const; + + //- Return access to mass + inline dimensionedScalar& mass(); + + //- Return moment of inertia + inline const dimensionedDiagTensor& momentOfInertia() const; + + //- Return access to moment of inertia + inline dimensionedDiagTensor& momentOfInertia(); + + //- Return equilibrium position of origin + inline const dimensionedVector& Xequilibrium() const; + + //- Return access to equilibrium position of origin + inline dimensionedVector& Xequilibrium(); + + + // Access to forces and moments + + //- Return force in inertial coord. sys. + inline const dimensionedVector& force() const; + + //- Return access to force in inertial coord. sys. + inline dimensionedVector& force(); + + //- Return moment in inertial coord. sys. + inline const dimensionedVector& moment() const; + + //- Return access to moment in inertial coord. sys. + inline dimensionedVector& moment(); + + //- Return force in relative coord. sys. + inline const dimensionedVector& forceRelative() const; + + //- Return access to force in relative coord. sys. + inline dimensionedVector& forceRelative(); + + //- Return moment in relative coord. sys. + inline const dimensionedVector& momentRelative() const; + + //- Return access to moment in relative coord. sys. + inline dimensionedVector& momentRelative(); + + //- Return total force in inertial coord. sys. + inline dimensionedVector forceTotal() const; + + //- Return total moment in inertial coord. sys. + inline dimensionedVector momentTotal() const; + + //- Relax the force (acceleration) using Aitkens or fixed relaxation + void relaxAcceleration + ( + const scalar minRelFactor, + const scalar maxRelFactor + ); + + + // Virtual interface for 6DOF motion state + + // Variables in relative coordinate system + + //- Return displacement in translated coordinate system + // relative to spring equilibrium + virtual const dimensionedVector& Xrel() const = 0; + + //- Return rotational velocity in relative coordinate system + virtual const dimensionedVector& omega() const = 0; + + + // Displacement and rotation in the absolute coordinate system + + //- Return position of origin in absolute coordinate system + virtual dimensionedVector X() const = 0; + + //- Return velocity of origin + virtual const dimensionedVector& U() const = 0; + + //- Return average velocity of origin for the previous time-step + virtual const dimensionedVector& Uaverage() const = 0; + + //- Return finite rotation + virtual const finiteRotation& rotation() const = 0; + + //- Return rotational vector of body + virtual vector rotVector() const = 0; + + //- Return rotation angle of body + virtual dimensionedScalar rotAngle() const = 0; + + + // Non-access control for setting state variables + + //- Set position of origin + virtual void setXrel(const vector& x) = 0; + + //- Set velocity of origin + virtual void setU(const vector& u) = 0; + + //- Set rotational angle in relative coordinate system + virtual void setRotation(const HamiltonRodriguezRot& rot) = 0; + + //- Set rotational velocity in relative coordinate system + virtual void setOmega(const vector& omega) = 0; + + //- Set referent coordinate system to apply constraints + virtual void setReferentRotation + ( + const HamiltonRodriguezRot& rot + ) = 0; + + + // Average motion per time-step + + //- Return average rotational vector of body + virtual vector rotVectorAverage() const = 0; + + //- Return average rotational velocity in relative coordinate + // system + virtual const dimensionedVector& omegaAverage() const = 0; + + //- Return average rotational velocity in absolute coordinate + // system + virtual const dimensionedVector& + omegaAverageAbsolute() const = 0; + + + // Rotations + + //- Return rotation tensor to relative coordinate system + virtual tensor toRelative() const = 0; + + //- Return rotation tensor to absolute coordinate system + virtual tensor toAbsolute() const = 0; + + //- Return transformation tensor between new and previous + // rotation + virtual const tensor& rotIncrementTensor() const = 0; + + //- Set ODE parameters from another ODE + virtual void setState(const sixDOFODE&) = 0; + + + // ODE parameters + + //- Return number of equations + virtual label nEqns() const = 0; + + //- Return access to coefficients + virtual scalarField& coeffs() = 0; + + //- Return reference to coefficients + virtual const scalarField& coeffs() const = 0; + + //- Evaluate derivatives + virtual void derivatives + ( + const scalar x, + const scalarField& y, + scalarField& dydx + ) const = 0; + + //- Evaluate Jacobian + virtual void jacobian + ( + const scalar x, + const scalarField& y, + scalarField& dfdx, + scalarSquareMatrix& dfdy + ) const = 0; + + //- Update ODE after the solution, advancing by delta + virtual void update(const scalar delta) = 0; + + + // Write control + + //- writeData member function required by regIOobject + virtual bool writeData(Ostream&) const; + + + // Ostream operator + + friend Ostream& operator<<(Ostream&, const sixDOFODE&); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "sixDOFODEI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/ODE/translationODE/translationODEI.H b/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H similarity index 51% rename from src/ODE/translationODE/translationODEI.H rename to src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H index 42c9ae94b..bb0d051d2 100644 --- a/src/ODE/translationODE/translationODEI.H +++ b/src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H @@ -21,74 +21,104 @@ License You should have received a copy of the GNU General Public License along with foam-extend. If not, see . -Class - translationODE - -Description - Ordinary differential equation for three degrees of freedom - solid body translation - -Author - Hrvoje Jasak - Dubravko Matijasevic - \*---------------------------------------------------------------------------*/ - // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -const Foam::dimensionedScalar& Foam::translationODE::mass() const +const Foam::dimensionedScalar& Foam::sixDOFODE::mass() const { return mass_; } -const Foam::dimensionedVector& Foam::translationODE::Xrel() const +Foam::dimensionedScalar& Foam::sixDOFODE::mass() { - return Xrel_; + return mass_; } -Foam::dimensionedVector Foam::translationODE::X() const +const Foam::dimensionedDiagTensor& Foam::sixDOFODE::momentOfInertia() const { - return xEquilibrium_ + Xrel(); + return momentOfInertia_; } -const Foam::dimensionedVector& Foam::translationODE::U() const +Foam::dimensionedDiagTensor& Foam::sixDOFODE::momentOfInertia() { - return U_; + return momentOfInertia_; } -const Foam::dimensionedVector& Foam::translationODE::Uold() const +const Foam::dimensionedVector& Foam::sixDOFODE::Xequilibrium() const { - return Uold_; + return Xequilibrium_; } -Foam::dimensionedVector Foam::translationODE::A() const +Foam::dimensionedVector& Foam::sixDOFODE::Xequilibrium() { - return A(Xrel(), U()); + return Xequilibrium_; } -Foam::dimensionedVector Foam::translationODE::Uaverage() const -{ - return 0.5*(U_ + Uold_); -} - - -const Foam::dimensionedVector& Foam::translationODE::force() const +const Foam::dimensionedVector& Foam::sixDOFODE::force() const { return force_; } -Foam::dimensionedVector& Foam::translationODE::force() +Foam::dimensionedVector& Foam::sixDOFODE::force() { return force_; } +const Foam::dimensionedVector& Foam::sixDOFODE::moment() const +{ + return moment_; +} + + +Foam::dimensionedVector& Foam::sixDOFODE::moment() +{ + return moment_; +} + + +const Foam::dimensionedVector& Foam::sixDOFODE::forceRelative() const +{ + return forceRelative_; +} + + +Foam::dimensionedVector& Foam::sixDOFODE::forceRelative() +{ + return forceRelative_; +} + + +const Foam::dimensionedVector& Foam::sixDOFODE::momentRelative() const +{ + return momentRelative_; +} + + +Foam::dimensionedVector& Foam::sixDOFODE::momentRelative() +{ + return momentRelative_; +} + + +Foam::dimensionedVector Foam::sixDOFODE::forceTotal() const +{ + return force_ + (this->toAbsolute() & forceRelative_); +} + + +Foam::dimensionedVector Foam::sixDOFODE::momentTotal() const +{ + return moment_ + (this->toAbsolute() & momentRelative_); +} + + // ************************************************************************* // diff --git a/src/ODE/translationODE/translationODE.C b/src/ODE/translationODE/translationODE.C deleted file mode 100644 index 099e87e66..000000000 --- a/src/ODE/translationODE/translationODE.C +++ /dev/null @@ -1,202 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend 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. - - foam-extend 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 foam-extend. If not, see . - -Class - translationODE - -Description - Ordinary differential equation for three degrees of freedom - solid body motion - -Author - Hrvoje Jasak - Dubravko Matijasevic - -\*---------------------------------------------------------------------------*/ - -#include "translationODE.H" -#include "objectRegistry.H" -#include "foamTime.H" - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -// namespace Foam -// { -// defineTypeNameAndDebug(translationODE, 0); -// } - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -void Foam::translationODE::setCoeffs() -{ - // Set ODE coefficients from position and rotation - - // Linear displacement in relative coordinate system - { - const vector& xVal = Xrel_.value(); - coeffs_[0] = xVal.x(); - coeffs_[1] = xVal.y(); - coeffs_[2] = xVal.z(); - } - - // Linear velocity in relative coordinate system - { - const vector& uVal = U_.value(); - coeffs_[3] = uVal.x(); - coeffs_[4] = uVal.y(); - coeffs_[5] = uVal.z(); - } -} - - -Foam::dimensionedVector Foam::translationODE::A -( - const dimensionedVector& xR, - const dimensionedVector& uR -) const -{ - return - ( - - (linSpringCoeffs_ & xR) // spring - - (linDampingCoeffs_ & uR) // damping - + force() - )/mass_; // gravity -} - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -// Construct from components -Foam::translationODE::translationODE -( - const IOobject& io -) -: - IOdictionary(io), - mass_(lookup("mass")), - xEquilibrium_(lookup("equilibriumPosition")), - linSpringCoeffs_(lookup("linearSpring")), - linDampingCoeffs_(lookup("linearDamping")), - Xrel_(lookup("Xrel")), - U_(lookup("U")), - Uold_(lookup("Uold")), - force_(lookup("force")), - coeffs_(6, 0.0) -{ - setCoeffs(); -} - - -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -Foam::translationODE::~translationODE() -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -void Foam::translationODE::derivatives -( - const scalar x, - const scalarField& y, - scalarField& dydx -) const -{ - // Set the derivatives for displacement - dydx[0] = y[3]; - dydx[1] = y[4]; - dydx[2] = y[5]; - - dimensionedVector curX("curX", dimLength, vector(y[0], y[1], y[2])); - dimensionedVector curU("curU", dimLength/dimTime, vector(y[3], y[4], y[5])); - - const vector& accel = A(curX, curU).value(); - - dydx[3] = accel.x(); - dydx[4] = accel.y(); - dydx[5] = accel.z(); -} - - -void Foam::translationODE::jacobian -( - const scalar x, - const scalarField& y, - scalarField& dfdx, - scalarSquareMatrix& dfdy -) const -{ - notImplemented("translationODE::jacobian(...) const"); -} - - -void Foam::translationODE::update(const scalar delta) -{ - // Update position - vector& Xval = Xrel_.value(); - - Xval.x() = coeffs_[0]; - Xval.y() = coeffs_[1]; - Xval.z() = coeffs_[2]; - - // Update velocity - Uold_ = U_; - - vector& Uval = U_.value(); - - Uval.x() = coeffs_[3]; - Uval.y() = coeffs_[4]; - Uval.z() = coeffs_[5]; -} - - -bool Foam::translationODE::writeData(Ostream& os) const -{ - os << *this; - return os.good(); -} - - -// * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * // - -Foam::Ostream& Foam::operator<<(Ostream& os, const translationODE& sds) -{ - os.writeKeyword("mass") << sds.mass_ << token::END_STATEMENT << endl; - os.writeKeyword("equilibriumPosition") << sds.xEquilibrium_ - << token::END_STATEMENT << endl; - os.writeKeyword("linearSpring") << sds.linSpringCoeffs_ - << token::END_STATEMENT << endl; - os.writeKeyword("linearDamping") << sds.linDampingCoeffs_ - << token::END_STATEMENT << endl; - - os.writeKeyword("Xrel") << sds.Xrel() << token::END_STATEMENT << endl; - os.writeKeyword("U") << sds.U() << token::END_STATEMENT << endl; - os.writeKeyword("Uold") << tab << sds.Uold() << token::END_STATEMENT << nl; - - os.writeKeyword("force") << sds.force() << token::END_STATEMENT << endl; - - return os; -} - - -// ************************************************************************* // diff --git a/src/ODE/translationODE/translationODE.H b/src/ODE/translationODE/translationODE.H deleted file mode 100644 index c66a479c9..000000000 --- a/src/ODE/translationODE/translationODE.H +++ /dev/null @@ -1,251 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | foam-extend: Open Source CFD - \\ / O peration | Version: 4.0 - \\ / A nd | Web: http://www.foam-extend.org - \\/ M anipulation | For copyright notice see file Copyright -------------------------------------------------------------------------------- -License - This file is part of foam-extend. - - foam-extend 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. - - foam-extend 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 foam-extend. If not, see . - -Class - translationODE - -Description - Ordinary differential equation for three degrees of freedom solid - body translation - -Author - Hrvoje Jasak - Dubravko Matijasevic - -SourceFiles - translationODEI.H - translationODE.C - -\*---------------------------------------------------------------------------*/ - -#ifndef translationODE_H -#define translationODE_H - -#include "ODE.H" -#include "IOdictionary.H" -#include "dimensionedTypes.H" -#include "dimensionedDiagTensor.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -namespace Foam -{ - -/*---------------------------------------------------------------------------*\ - Class translationODE Declaration -\*---------------------------------------------------------------------------*/ - -class translationODE -: - public IOdictionary, - public ODE -{ - // Private data - - // Body data - - //- Mass - dimensionedScalar mass_; - - - // Platform variables - - //- Spring equilibrium position for translation - const dimensionedVector xEquilibrium_; - - //- Linear spring coeffs - const dimensionedDiagTensor linSpringCoeffs_; - - //- Linear damping coeffs - const dimensionedDiagTensor linDampingCoeffs_; - - - // Body position and rotation variables - - //- Displacement relative to spring equilibrium - dimensionedVector Xrel_; - - //- Velocity of mass centroid - dimensionedVector U_; - - //- Velocity of mass centroid at previous time-step - dimensionedVector Uold_; - - - // External forces - - //- Force driving the motion - dimensionedVector force_; - - - //- ODE coefficients - scalarField coeffs_; - - - // Private Member Functions - - //- Disallow default bitwise copy construct - translationODE(const translationODE&); - - //- Disallow default bitwise assignment - void operator=(const translationODE&); - - - //- Set ODE coefficients from position and rotation - inline void setCoeffs(); - - - // Variables in relative coordinate system (solved for) - - //- Return acceleration in relative coordinate system - // given current values of relative displacement and velocity - dimensionedVector A - ( - const dimensionedVector& xR, - const dimensionedVector& uR - ) const; - - -public: - -// //- Runtime type information -// TypeName("translationODE"); - - - // Constructors - - //- Construct from dictionary - translationODE(const IOobject& io); - - - // Destructor - - virtual ~translationODE(); - - - // Member Functions - - //- Return mass - inline const dimensionedScalar& mass() const; - - - // Variables in relative coordinate system (solved for) - - //- Return displacement in relative coordinate system - inline const dimensionedVector& Xrel() const; - - - // Displacement and rotation in the absolute coordinate system - - //- Return position of origin in absolute coordinate system - inline dimensionedVector X() const; - - //- Return velocity of origin - inline const dimensionedVector& U() const; - - //- Return velocity of origin for the previous time-step - inline const dimensionedVector& Uold() const; - - //- Return acceleration of origin - inline dimensionedVector A() const; - - - // Average motion per time-step - - //- Return average velocity of origin - inline dimensionedVector Uaverage() const; - - - // Force - - //- Return force - inline const dimensionedVector& force() const; - - //- Return access to force - inline dimensionedVector& force(); - - - // ODE parameters - - //- Return number of equations - virtual label nEqns() const - { - return 6; - } - - //- Return reference to interpolation coefficients - virtual scalarField& coeffs() - { - return coeffs_; - } - - //- Return reference to interpolation coefficients - virtual const scalarField& coeffs() const - { - return coeffs_; - } - - //- Return derivatives - virtual void derivatives - ( - const scalar x, - const scalarField& y, - scalarField& dydx - ) const; - - //- Return Jacobian - virtual void jacobian - ( - const scalar x, - const scalarField& y, - scalarField& dfdx, - scalarSquareMatrix& dfdy - ) const; - - //- Update ODE after the solution, advancing by delta - virtual void update(const scalar delta); - - - //- WriteData member function required by regIOobject - bool writeData(Ostream&) const; - - - // Ostream operator - - friend Ostream& operator<<(Ostream&, const translationODE&); -}; - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -} // End namespace Foam - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#include "translationODEI.H" - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* //