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

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2010-2011 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "tabulatedAxialAngularSpring.H"
  27 #include "addToRunTimeSelectionTable.H"
  28 #include "sixDoFRigidBodyMotion.H"
  29 #include "transform.H"
  30 #include "unitConversion.H"
  31
  32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36 namespace sixDoFRigidBodyMotionRestraints
  37 {
  38     defineTypeNameAndDebug(tabulatedAxialAngularSpring, 0);
  39
  40     addToRunTimeSelectionTable
  41     (
  42         sixDoFRigidBodyMotionRestraint,
  43         tabulatedAxialAngularSpring,
  44         dictionary
  45     );
  46 }
  47 }
  48
  49
  50 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  51
  52 Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring::
  53 tabulatedAxialAngularSpring
  54 (
  55     const dictionary& sDoFRBMRDict
  56 )
  57 :
  58     sixDoFRigidBodyMotionRestraint(sDoFRBMRDict),
  59     refQ_(),
  60     axis_(),
  61     moment_(),
  62     convertToDegrees_(),
  63     damping_()
  64 {
  65     read(sDoFRBMRDict);
  66 }
  67
  68
  69 // * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
  70
  71 Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring::
  72 ~tabulatedAxialAngularSpring()
  73 {}
  74
  75
  76 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
  77
  78 void
  79 Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring::restrain
  80 (
  81     const sixDoFRigidBodyMotion& motion,
  82     vector& restraintPosition,
  83     vector& restraintForce,
  84     vector& restraintMoment
  85 ) const
  86 {
  87     vector refDir = rotationTensor(vector(1, 0 ,0), axis_) & vector(0, 1, 0);
  88
  89     vector oldDir = refQ_ & refDir;
  90     vector newDir = motion.orientation() & refDir;
  91
  92     if (mag(oldDir & axis_) > 0.95 || mag(newDir & axis_) > 0.95)
  93     {
  94         // Directions getting close to the axis, change reference
  95
  96         refDir = rotationTensor(vector(1, 0 ,0), axis_) & vector(0, 0, 1);
  97
  98         oldDir = refQ_ & refDir;
  99         newDir = motion.orientation() & refDir;
 100     }
 101
 102     // Removing any axis component from oldDir and newDir and normalising
 103     oldDir -= (axis_ & oldDir)*axis_;
 104     oldDir /= (mag(oldDir) + VSMALL);
 105
 106     newDir -= (axis_ & newDir)*axis_;
 107     newDir /= (mag(newDir) + VSMALL);
 108
 109     scalar theta = mag(acos(min(oldDir & newDir, 1.0)));
 110
 111     // Determining the sign of theta by comparing the cross product of
 112     // the direction vectors with the axis
 113     theta *= sign((oldDir ^ newDir) & axis_);
 114
 115     scalar moment;
 116
 117     if (convertToDegrees_)
 118     {
 119         moment = moment_(radToDeg(theta));
 120     }
 121     else
 122     {
 123         moment = moment_(theta);
 124     }
 125
 126     // Damping of along axis angular velocity only
 127     restraintMoment = moment*axis_ - damping_*(motion.omega() & axis_)*axis_;
 128
 129     restraintForce = vector::zero;
 130
 131     // Not needed to be altered as restraintForce is zero, but set to
 132     // centreOfMass to be sure of no spurious moment
 133     restraintPosition = motion.centreOfMass();
 134
 135     if (motion.report())
 136     {
 137         Info<< " angle " << theta
 138             << " force " << restraintForce
 139             << " moment " << restraintMoment
 140             << endl;
 141     }
 142 }
 143
 144
 145 bool Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring::read
 146 (
 147     const dictionary& sDoFRBMRDict
 148 )
 149 {
 150     sixDoFRigidBodyMotionRestraint::read(sDoFRBMRDict);
 151
 152     refQ_ = sDoFRBMRCoeffs_.lookupOrDefault<tensor>("referenceOrientation", I);
 153
 154     if (mag(mag(refQ_) - sqrt(3.0)) > 1e-9)
 155     {
 156         FatalErrorIn
 157         (
 158             "Foam::sixDoFRigidBodyMotionRestraints::"
 159             "tabulatedAxialAngularSpring::read"
 160             "("
 161                 "const dictionary& sDoFRBMRDict"
 162             ")"
 163         )
 164             << "referenceOrientation " << refQ_ << " is not a rotation tensor. "
 165             << "mag(referenceOrientation) - sqrt(3) = "
 166             << mag(refQ_) - sqrt(3.0) << nl
 167             << exit(FatalError);
 168     }
 169
 170     axis_ = sDoFRBMRCoeffs_.lookup("axis");
 171
 172     scalar magAxis(mag(axis_));
 173
 174     if (magAxis > VSMALL)
 175     {
 176         axis_ /= magAxis;
 177     }
 178     else
 179     {
 180         FatalErrorIn
 181         (
 182             "Foam::sixDoFRigidBodyMotionRestraints::"
 183             "tabulatedAxialAngularSpring::read"
 184             "("
 185                 "const dictionary& sDoFRBMCDict"
 186             ")"
 187         )
 188             << "axis has zero length"
 189             << abort(FatalError);
 190     }
 191
 192     moment_ = interpolationTable<scalar>(sDoFRBMRCoeffs_);
 193
 194     const word angleFormat = sDoFRBMRCoeffs_.lookup("angleFormat");
 195
 196     if (angleFormat == "degrees" || angleFormat == "degree")
 197     {
 198         convertToDegrees_ = true;
 199     }
 200     else if (angleFormat == "radians" || angleFormat == "radian")
 201     {
 202         convertToDegrees_ = false;
 203     }
 204     else
 205     {
 206         FatalErrorIn
 207         (
 208             "Foam::sixDoFRigidBodyMotionRestraints::"
 209             "tabulatedAxialAngularSpring::read"
 210             "("
 211                 "const dictionary&"
 212             ")"
 213         )
 214             << "angleFormat must be degree, degrees, radian or radians"
 215             << abort(FatalError);
 216     }
 217
 218     sDoFRBMRCoeffs_.lookup("damping") >> damping_;
 219
 220     return true;
 221 }
 222
 223
 224 void Foam::sixDoFRigidBodyMotionRestraints::tabulatedAxialAngularSpring::write
 225 (
 226     Ostream& os
 227 ) const
 228 {
 229     os.writeKeyword("referenceOrientation")
 230         << refQ_ << token::END_STATEMENT << nl;
 231
 232     os.writeKeyword("axis")
 233         << axis_ << token::END_STATEMENT << nl;
 234
 235     moment_.write(os);
 236
 237     os.writeKeyword("angleFormat");
 238
 239     if (convertToDegrees_)
 240     {
 241         os  << "degrees" << token::END_STATEMENT << nl;
 242     }
 243     else
 244     {
 245         os  << "radians" << token::END_STATEMENT << nl;
 246     }
 247
 248     os.writeKeyword("damping")
 249         << damping_ << token::END_STATEMENT << nl;
 250 }
 251
 252
 253 // ************************************************************************* //