ENH: sixDoFRigidBodyMotion. Adding linearSphericalAngularSpring restraint.

Moving maximum iteration FatalErrorIn check after constraint while loop.

Adding currentOrientation function to give access to Q, but safely and
with a more descriptive name for the function.

src/postProcessing/functionObjects/forces/Make/files

@@ -16,6 +16,7 @@ sDoFRBMR = $(sDoFRBM)/sixDoFRigidBodyMotionRestraint
 $(sDoFRBMR)/sixDoFRigidBodyMotionRestraint/sixDoFRigidBodyMotionRestraint.C
 $(sDoFRBMR)/sixDoFRigidBodyMotionRestraint/newSixDoFRigidBodyMotionRestraint.C
 $(sDoFRBMR)/linearSpring/linearSpring.C
+$(sDoFRBMR)/linearSphericalAngularSpring/linearSphericalAngularSpring.C
 
 sDoFRBMC = $(sDoFRBM)/sixDoFRigidBodyMotionConstraint
 

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C

@@ -45,7 +45,7 @@ void Foam::sixDoFRigidBodyMotion::applyRestraints()
 
             restraints_[rI].restrain(*this, rP, rF, rM);
 
-            Info<< "Restraint " << rI << " force " << rF << endl;
+            // Info<< "Restraint " << rI << " force " << rF << endl;
 
             a() += rF/mass_;
 
@@ -71,26 +71,11 @@ void Foam::sixDoFRigidBodyMotion::applyConstraints(scalar deltaT)
 
         do
         {
-            iter++;
-
-            if (iter > maxConstraintIters_)
-            {
-                FatalErrorIn
-                (
-                    "Foam::sixDoFRigidBodyMotion::applyConstraints"
-                    "(scalar deltaT)"
-                )
-                    << nl << "Maximum number of constraint iterations ("
-                    << maxConstraintIters_ << ") exceeded." << nl
-                    << exit(FatalError);
-
-            }
-
             allConverged = true;
 
             forAll(constraints_, cI)
             {
-                Info<< nl << "Constraint " << cI << endl;
+                // Info<< nl << "Constraint " << cI << endl;
 
                 // constraint position
                 point cP = vector::zero;
@@ -121,12 +106,27 @@ void Foam::sixDoFRigidBodyMotion::applyConstraints(scalar deltaT)
                 cMA += cM + ((cP - centreOfMass()) ^ cF);
             }
 
-        } while(!allConverged);
+        } while(!allConverged && ++iter < maxConstraintIters_);
 
-        Info<< "Constraints converged in " << iter << " iterations" << nl
+        if (iter >= maxConstraintIters_)
-            << "Constraint force: " << cFA << nl
+        {
-            << "Constraint moment: " << cMA
+            FatalErrorIn
+            (
+                "Foam::sixDoFRigidBodyMotion::applyConstraints"
+                "(scalar deltaT)"
+            )
+                << nl << "Maximum number of sixDoFRigidBodyMotion constraint "
+                << "iterations (" << maxConstraintIters_ << ") exceeded." << nl
+                << exit(FatalError);
+        }
+        else
+        {
+            Info<< "sixDoFRigidBodyMotion constraints converged in "
+                << iter << " iterations" << nl
+            // << "Constraint force: " << cFA << nl
+            // << "Constraint moment: " << cMA
                 << endl;
+        }
 
         // Add the constrain forces and moments to the motion state variables
         a() += cFA/mass_;

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.H

@@ -272,6 +272,10 @@ public:
         //  motion state
         inline point currentPosition(const point& refPt) const;
 
+        //- Transform the given reference state direction by the current
+        //  motion state
+        inline point currentOrientation(const vector& refDir) const;
+
         //- Predict the position of the supplied point after deltaT
         //  given the current motion state and the additional supplied
         //  force and moment

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedLine/fixedLine.C

@@ -89,16 +89,16 @@ bool Foam::sixDoFRigidBodyMotionConstraints::fixedLine::constrain
 
     constraintPosition = motion.currentPosition(refPt_);
 
-    Info<< "current position " << constraintPosition << nl
+    // Info<< "current position " << constraintPosition << nl
-        << "next predictedPosition " << predictedPosition
+    //     << "next predictedPosition " << predictedPosition
-        << endl;
+    //     << endl;
 
     // Vector from reference point to predicted point
     vector rC = predictedPosition - refPt_;
 
     vector error = rC - ((rC) & dir_)*dir_;
 
-    Info<< "error " << error << endl;
+    // Info<< "error " << error << endl;
 
     constraintForceIncrement =
         -relaxationFactor_*error*motion.mass()/sqr(deltaT);

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionConstraint/fixedPlane/fixedPlane.C

@@ -92,13 +92,13 @@ bool Foam::sixDoFRigidBodyMotionConstraints::fixedPlane::constrain
 
     constraintPosition = motion.currentPosition(refPt);
 
-    Info<< "current position " << constraintPosition << nl
+    // Info<< "current position " << constraintPosition << nl
-        << "next predictedPosition " << predictedPosition
+    //     << "next predictedPosition " << predictedPosition
-        << endl;
+    //     << endl;
 
     vector error = ((predictedPosition - refPt) & n)*n;
 
-    Info<< "error " << error << endl;
+    // Info<< "error " << error << endl;
 
     constraintForceIncrement =
         -relaxationFactor_*error*motion.mass()/sqr(deltaT);

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionI.H

@@ -223,6 +223,15 @@ inline Foam::point Foam::sixDoFRigidBodyMotion::currentPosition
 }
 
 
+inline Foam::point Foam::sixDoFRigidBodyMotion::currentOrientation
+(
+    const vector& refDir
+) const
+{
+    return (Q() & refDir);
+}
+
+
 inline Foam::vector Foam::sixDoFRigidBodyMotion::omega() const
 {
     return  Q() & (inv(momentOfInertia_) & pi());

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSphericalAngularSpring/linearSphericalAngularSpring.C

@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 "linearSphericalAngularSpring.H"
+#include "addToRunTimeSelectionTable.H"
+#include "sixDoFRigidBodyMotion.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace sixDoFRigidBodyMotionRestraints
+{
+    defineTypeNameAndDebug(linearSphericalAngularSpring, 0);
+    addToRunTimeSelectionTable
+    (
+        sixDoFRigidBodyMotionRestraint,
+        linearSphericalAngularSpring,
+        dictionary
+    );
+};
+};
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sixDoFRigidBodyMotionRestraints::linearSphericalAngularSpring::
+linearSphericalAngularSpring
+(
+    const dictionary& sDoFRBMRDict
+)
+:
+    sixDoFRigidBodyMotionRestraint(sDoFRBMRDict),
+    refDir_(),
+    stiffness_(),
+    damping_()
+{
+    read(sDoFRBMRDict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
+
+Foam::sixDoFRigidBodyMotionRestraints::linearSphericalAngularSpring::
+~linearSphericalAngularSpring()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void
+Foam::sixDoFRigidBodyMotionRestraints::linearSphericalAngularSpring::restrain
+(
+    const sixDoFRigidBodyMotion& motion,
+    vector& restraintPosition,
+    vector& restraintForce,
+    vector& restraintMoment
+) const
+{
+    vector curDir = motion.currentOrientation(refDir_);
+
+    vector a = refDir_ ^ curDir;
+
+    scalar magA = mag(a);
+
+    scalar theta = 0;
+
+    if (magA > VSMALL)
+    {
+        a /= magA;
+
+        theta = acos(curDir & refDir_);
+    }
+    else
+    {
+        a = vector::zero;
+    }
+
+    restraintMoment = -stiffness_*theta*a - damping_*motion.omega();
+
+    restraintForce = vector::zero;
+
+    // Not needed to be altered as restraintForce is zero, but set to
+    // centre of mass to be sure of no spurious moment
+    restraintPosition = motion.centreOfMass();
+}
+
+
+bool Foam::sixDoFRigidBodyMotionRestraints::linearSphericalAngularSpring ::read
+(
+    const dictionary& sDoFRBMRDict
+)
+{
+    sixDoFRigidBodyMotionRestraint::read(sDoFRBMRDict);
+
+    sDoFRBMRCoeffs_.lookup("referenceDirection") >> refDir_;
+
+    sDoFRBMRCoeffs_.lookup("stiffness") >> stiffness_;
+
+    sDoFRBMRCoeffs_.lookup("damping") >> damping_;
+
+    return true;
+}
+
+// ************************************************************************* //

src/postProcessing/functionObjects/forces/pointPatchFields/derived/sixDoFRigidBodyMotion/sixDoFRigidBodyMotionRestraint/linearSphericalAngularSpring/linearSphericalAngularSpring.H

@@ -0,0 +1,125 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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::sixDoFRigidBodyMotionRestraints::linearSphericalAngularSpring
+
+Description
+    sixDoFRigidBodyMotionRestraints model.  Linear spherical angular spring.
+
+SourceFiles
+    linearSphericalAngularSpring.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef linearSphericalAngularSpring_H
+#define linearSphericalAngularSpring_H
+
+#include "sixDoFRigidBodyMotionRestraint.H"
+#include "point.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace sixDoFRigidBodyMotionRestraints
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class linearSphericalAngularSpring Declaration
+\*---------------------------------------------------------------------------*/
+
+class linearSphericalAngularSpring
+:
+    public sixDoFRigidBodyMotionRestraint
+{
+    // Private data
+
+        // Reference direction there is no spring reaction torque
+        vector refDir_;
+
+        //- Spring stiffness coefficient (Nm/rad)
+        scalar stiffness_;
+
+        //- Damping coefficient (Nms/rad)
+        scalar damping_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("linearSphericalAngularSpring");
+
+
+    // Constructors
+
+        //- Construct from components
+        linearSphericalAngularSpring
+        (
+            const dictionary& sDoFRBMRDict
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<sixDoFRigidBodyMotionRestraint> clone() const
+        {
+            return autoPtr<sixDoFRigidBodyMotionRestraint>
+            (
+                new linearSphericalAngularSpring(*this)
+            );
+        }
+
+
+    // Destructor
+
+        virtual ~linearSphericalAngularSpring();
+
+
+    // Member Functions
+
+        //- Calculate the restraint position, force and moment.
+        //  Global reference frame vectors.
+        virtual void restrain
+        (
+            const sixDoFRigidBodyMotion& motion,
+            vector& restraintPosition,
+            vector& restraintForce,
+            vector& restraintMoment
+        ) const;
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& sDoFRBMRCoeff);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace solidBodyMotionFunctions
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //