diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculateExternalForce.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculateExternalForce.C index 2def748e1d..20f2776e04 100644 --- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculateExternalForce.C +++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculateExternalForce.C @@ -32,9 +32,13 @@ void Foam::moleculeCloud::calculateExternalForce() { iterator mol(this->begin()); + Info<< "Warning! Includes dissipation term!" << endl; + for (mol = this->begin(); mol != this->end(); ++mol) { mol().A() += gravity_; + + mol().A() += -1.0 * mol().U() /mol().mass(); } } diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculateTetherForce.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculateTetherForce.C index da2dd57121..3347a4e041 100644 --- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculateTetherForce.C +++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudCalculateTetherForce.C @@ -34,8 +34,6 @@ void Foam::moleculeCloud::calculateTetherForce() vector rIT; - scalar rITMag; - vector fIT; for (mol = this->begin(); mol != this->end(); ++mol) @@ -44,12 +42,10 @@ void Foam::moleculeCloud::calculateTetherForce() { rIT = mol().position() - mol().tetherPosition(); - rITMag = mag(rIT); - - fIT = (rIT/rITMag) * tetherPotentials_.force + fIT = tetherPotentials_.force ( mol().id(), - rITMag + rIT ); mol().A() += fIT/(mol().mass()); @@ -57,7 +53,7 @@ void Foam::moleculeCloud::calculateTetherForce() mol().potentialEnergy() += tetherPotentials_.energy ( mol().id(), - rITMag + rIT ); mol().rf() += rIT*fIT; diff --git a/src/lagrangian/molecularDynamics/potential/Make/files b/src/lagrangian/molecularDynamics/potential/Make/files index 43cf0bb737..1e2b01c23a 100644 --- a/src/lagrangian/molecularDynamics/potential/Make/files +++ b/src/lagrangian/molecularDynamics/potential/Make/files @@ -9,6 +9,7 @@ $(pairPotential)/basic/newPairPotential.C $(pairPotential)/derived/lennardJones/lennardJones.C $(pairPotential)/derived/maitlandSmith/maitlandSmith.C $(pairPotential)/derived/azizChen/azizChen.C +$(pairPotential)/derived/exponentialRepulsion/exponentialRepulsion.C energyScalingFunction = energyScalingFunction @@ -30,6 +31,7 @@ $(tetherPotential)/basic/newTetherPotential.C $(tetherPotential)/derived/harmonicSpring/harmonicSpring.C $(tetherPotential)/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.C +$(tetherPotential)/derived/pitchForkRing/pitchForkRing.C LIB = $(FOAM_LIBBIN)/libpotential diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/exponentialRepulsion/exponentialRepulsion.C b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/exponentialRepulsion/exponentialRepulsion.C new file mode 100644 index 0000000000..3034834ccf --- /dev/null +++ b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/exponentialRepulsion/exponentialRepulsion.C @@ -0,0 +1,93 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ 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 "exponentialRepulsion.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace pairPotentials +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(exponentialRepulsion, 0); + +addToRunTimeSelectionTable +( + pairPotential, + exponentialRepulsion, + dictionary +); + +// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +exponentialRepulsion::exponentialRepulsion +( + const word& name, + const dictionary& exponentialRepulsion +) +: + pairPotential(name, exponentialRepulsion), + exponentialRepulsionCoeffs_(exponentialRepulsion.subDict(typeName + "Coeffs")), + rm_(readScalar(exponentialRepulsionCoeffs_.lookup("rm"))), + epsilon_(readScalar(exponentialRepulsionCoeffs_.lookup("epsilon"))) +{ + setLookupTables(); +} + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +scalar exponentialRepulsion::unscaledEnergy(const scalar r) const +{ + return epsilon_ * exp(-r/rm_); +} + + +bool exponentialRepulsion::read(const dictionary& exponentialRepulsion) +{ + pairPotential::read(exponentialRepulsion); + + exponentialRepulsionCoeffs_ = exponentialRepulsion.subDict(typeName + "Coeffs"); + + exponentialRepulsionCoeffs_.lookup("rm") >> rm_; + exponentialRepulsionCoeffs_.lookup("epsilon") >> epsilon_; + + return true; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace pairPotentials +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/lagrangian/molecularDynamics/potential/pairPotential/derived/exponentialRepulsion/exponentialRepulsion.H b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/exponentialRepulsion/exponentialRepulsion.H new file mode 100644 index 0000000000..7026710427 --- /dev/null +++ b/src/lagrangian/molecularDynamics/potential/pairPotential/derived/exponentialRepulsion/exponentialRepulsion.H @@ -0,0 +1,103 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ 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 + +Class + Foam::pairPotentials::exponentialRepulsion + +Description + +SourceFiles + exponentialRepulsion.C + +\*---------------------------------------------------------------------------*/ + +#ifndef exponentialRepulsion_H +#define exponentialRepulsion_H + +#include "pairPotential.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +namespace pairPotentials +{ + +/*---------------------------------------------------------------------------*\ + Class exponentialRepulsion Declaration +\*---------------------------------------------------------------------------*/ + +class exponentialRepulsion +: + public pairPotential +{ + // Private data + + dictionary exponentialRepulsionCoeffs_; + + scalar rm_; + scalar epsilon_; + +public: + + //- Runtime type information + TypeName("exponentialRepulsion"); + + + // Constructors + + //- Construct from components + exponentialRepulsion + ( + const word& name, + const dictionary& pairPotentialProperties + ); + + + // Destructor + + ~exponentialRepulsion() + {} + + + // Member Functions + + scalar unscaledEnergy(const scalar r) const; + + //- Read transportProperties dictionary + bool read(const dictionary& pairPotentialProperties); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace pairPotentials +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.H index e7ffb4e7a3..10ef201a14 100644 --- a/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.H +++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/basic/tetherPotential.H @@ -40,6 +40,7 @@ SourceFiles #include "typeInfo.H" #include "runTimeSelectionTables.H" #include "autoPtr.H" +#include "vector.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -119,9 +120,9 @@ public: // Member Functions - virtual scalar energy (const scalar r) const = 0; + virtual scalar energy (const vector r) const = 0; - virtual scalar force (const scalar r) const = 0; + virtual vector force (const vector r) const = 0; const dictionary& tetherPotentialProperties() const { diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.C b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.C index e21f077a1b..c26fa022ee 100644 --- a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.C +++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.C @@ -65,13 +65,13 @@ harmonicSpring::harmonicSpring // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // -scalar harmonicSpring::energy(const scalar r) const +scalar harmonicSpring::energy(const vector r) const { - return 0.5*springConstant_*r*r; + return 0.5*springConstant_*magSqr(r); } -scalar harmonicSpring::force(const scalar r) const +vector harmonicSpring::force(const vector r) const { return -springConstant_*r; } diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.H index 0c5057285c..774e4e4b72 100644 --- a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.H +++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/harmonicSpring/harmonicSpring.H @@ -84,9 +84,9 @@ public: // Member Functions - scalar energy(const scalar r) const; + scalar energy(const vector r) const; - scalar force(const scalar r) const; + vector force(const vector r) const; //- Read transportProperties dictionary bool read(const dictionary& tetherPotentialProperties); diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.C b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.C new file mode 100644 index 0000000000..87b9c7b5bd --- /dev/null +++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.C @@ -0,0 +1,117 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ 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 "pitchForkRing.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace tetherPotentials +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(pitchForkRing, 0); + +addToRunTimeSelectionTable +( + tetherPotential, + pitchForkRing, + dictionary +); + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +pitchForkRing::pitchForkRing +( + const word& name, + const dictionary& tetherPotentialProperties +) +: + tetherPotential(name, tetherPotentialProperties), + pitchForkRingCoeffs_ + ( + tetherPotentialProperties.subDict(typeName + "Coeffs") + ), + mu_(readScalar(pitchForkRingCoeffs_.lookup("mu"))), + alpha_(readScalar(pitchForkRingCoeffs_.lookup("alpha"))), + rOrbit_(readScalar(pitchForkRingCoeffs_.lookup("rOrbit"))) +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +scalar pitchForkRing::energy(const vector r) const +{ + scalar p = sqrt(r.x()*r.x() + r.y()*r.y()); + + scalar pMinusRSqr = (p - rOrbit_)*(p - rOrbit_); + + return -0.5 * mu_ * pMinusRSqr + + 0.25 * pMinusRSqr * pMinusRSqr + + 0.5 * alpha_ * r.z() * r.z(); +} + + +vector pitchForkRing::force(const vector r) const +{ + scalar p = sqrt(r.x()*r.x() + r.y()*r.y()); + + scalar pMinusR = (p - rOrbit_); + + return vector + ( + (mu_ - pMinusR * pMinusR) * pMinusR * r.x()/p, + (mu_ - pMinusR * pMinusR) * pMinusR * r.y()/p, + -alpha_ * r.z() + ); +} + + +bool pitchForkRing::read(const dictionary& tetherPotentialProperties) +{ + tetherPotential::read(tetherPotentialProperties); + + pitchForkRingCoeffs_ = + tetherPotentialProperties.subDict(typeName + "Coeffs"); + + pitchForkRingCoeffs_.lookup("mu") >> mu_; + pitchForkRingCoeffs_.lookup("alpha") >> alpha_; + pitchForkRingCoeffs_.lookup("rOrbit") >> rOrbit_; + + return true; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace tetherPotentials +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.H new file mode 100644 index 0000000000..25bb2d5e21 --- /dev/null +++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.H @@ -0,0 +1,107 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd. + \\/ 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 + +Class + Foam::tetherPotentials::pitchForkRing + +Description + + +SourceFiles + pitchForkRing.C + +\*---------------------------------------------------------------------------*/ + +#ifndef pitchForkRing_H +#define pitchForkRing_H + +#include "tetherPotential.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace tetherPotentials +{ + +/*---------------------------------------------------------------------------*\ + Class pitchForkRing Declaration +\*---------------------------------------------------------------------------*/ + +class pitchForkRing +: + public tetherPotential +{ + // Private data + + dictionary pitchForkRingCoeffs_; + + scalar mu_; + scalar alpha_; + scalar rOrbit_; + + +public: + + //- Runtime type information + TypeName("pitchForkRing"); + + + // Constructors + + //- Construct from components + pitchForkRing + ( + const word& name, + const dictionary& tetherPotentialProperties + ); + + + // Destructor + + ~pitchForkRing() + {} + + + // Member Functions + + scalar energy(const vector r) const; + + vector force(const vector r) const; + + //- Read transportProperties dictionary + bool read(const dictionary& tetherPotentialProperties); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace tetherPotentials +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.C b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.C index 05cf2ecda0..7a08b3ea39 100644 --- a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.C +++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.C @@ -73,30 +73,33 @@ restrainedHarmonicSpring::restrainedHarmonicSpring // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // -scalar restrainedHarmonicSpring::energy(const scalar r) const +scalar restrainedHarmonicSpring::energy(const vector r) const { - if (r < rR_) + scalar magR = mag(r); + + if (magR < rR_) { - return 0.5 * springConstant_ * r * r; + return 0.5 * springConstant_ * magSqr(r); } else { return 0.5 * springConstant_ * rR_ * rR_ - + springConstant_ * rR_ * (r - rR_); + + springConstant_ * rR_ * (magR - rR_); } } -scalar restrainedHarmonicSpring::force(const scalar r) const +vector restrainedHarmonicSpring::force(const vector r) const { - if (r < rR_) + scalar magR = mag(r); + + if (magR < rR_) { return -springConstant_ * r; } else { - return -springConstant_ * rR_; + return -springConstant_ * rR_ * r / magR; } - } bool restrainedHarmonicSpring::read(const dictionary& tetherPotentialProperties) @@ -108,7 +111,7 @@ bool restrainedHarmonicSpring::read(const dictionary& tetherPotentialProperties) restrainedHarmonicSpringCoeffs_.lookup("springConstant") >> springConstant_; restrainedHarmonicSpringCoeffs_.lookup("rR") >> rR_; - + return true; } diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.H index 185ba7f1ba..49880235ad 100644 --- a/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.H +++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/derived/restrainedHarmonicSpring/restrainedHarmonicSpring.H @@ -86,9 +86,9 @@ public: // Member Functions - scalar energy(const scalar r) const; + scalar energy(const vector r) const; - scalar force(const scalar r) const; + vector force(const vector r) const; //- Read transportProperties dictionary bool read(const dictionary& tetherPotentialProperties); diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.C b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.C index 1d5edf0c42..ee0d417b57 100644 --- a/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.C +++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.C @@ -136,27 +136,23 @@ const Foam::tetherPotential& Foam::tetherPotentialList::tetherPotentialFunction } -Foam::scalar Foam::tetherPotentialList::force +Foam::vector Foam::tetherPotentialList::force ( const label a, - const scalar rITMag + const vector rIT ) const { - scalar f = (*this)[tetherPotentialIndex (a)].force(rITMag); - - return f; + return (*this)[tetherPotentialIndex (a)].force(rIT); } Foam::scalar Foam::tetherPotentialList::energy ( const label a, - const scalar rITMag + const vector rIT ) const { - scalar e = (*this)[tetherPotentialIndex (a)].energy(rITMag); - - return e; + return (*this)[tetherPotentialIndex (a)].energy(rIT); } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.H b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.H index fb3fcc3758..3917d9329d 100644 --- a/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.H +++ b/src/lagrangian/molecularDynamics/potential/tetherPotential/tetherPotentialList/tetherPotentialList.H @@ -113,16 +113,16 @@ public: const label a ) const; - scalar force + vector force ( const label a, - const scalar rIJMag + const vector rIT ) const; scalar energy ( const label a, - const scalar rIJMag + const vector rIT ) const; };