linearAxialAngularSpring.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
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 
28 #include "sixDoFRigidBodyMotion.H"
29 #include "transform.H"
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35 namespace sixDoFRigidBodyMotionRestraints
36 {
38 
40  (
44  );
45 }
46 }
47 
48 
49 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
50 
53 (
54  const word& name,
55  const dictionary& sDoFRBMRDict
56 )
57 :
58  sixDoFRigidBodyMotionRestraint(name, sDoFRBMRDict),
59  refQ_(),
60  axis_(),
61  stiffness_(),
62  damping_()
63 {
64  read(sDoFRBMRDict);
65 }
66 
67 
68 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
69 
72 {}
73 
74 
75 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
76 
77 void
79 (
80  const sixDoFRigidBodyMotion& motion,
81  vector& restraintPosition,
82  vector& restraintForce,
83  vector& restraintMoment
84 ) const
85 {
86  vector refDir = rotationTensor(vector(1, 0 ,0), axis_) & vector(0, 1, 0);
87 
88  vector oldDir = refQ_ & refDir;
89  vector newDir = motion.orientation() & refDir;
90 
91  if (mag(oldDir & axis_) > 0.95 || mag(newDir & axis_) > 0.95)
92  {
93  // Directions getting close to the axis, change reference
94 
95  refDir = rotationTensor(vector(1, 0 ,0), axis_) & vector(0, 0, 1);
96  oldDir = refQ_ & refDir;
97  newDir = motion.orientation() & refDir;
98  }
99 
100  // Removing any axis component from oldDir and newDir and normalising
101  oldDir -= (axis_ & oldDir)*axis_;
102  oldDir /= (mag(oldDir) + VSMALL);
103 
104  newDir -= (axis_ & newDir)*axis_;
105  newDir /= (mag(newDir) + VSMALL);
106 
107  scalar theta = mag(acos(min(oldDir & newDir, 1.0)));
108 
109  // Temporary axis with sign information.
110  vector a = (oldDir ^ newDir);
111 
112  // Remove any component that is not along axis that may creep in
113  a = (a & axis_)*axis_;
114 
115  scalar magA = mag(a);
116 
117  if (magA > VSMALL)
118  {
119  a /= magA;
120  }
121  else
122  {
123  a = Zero;
124  }
125 
126  // Damping of along axis angular velocity only
127  restraintMoment = -stiffness_*theta*a - damping_*(motion.omega() & a)*a;
128 
129  restraintForce = Zero;
130 
131  // Not needed to be altered as restraintForce is zero, but set to
132  // centreOfRotation to be sure of no spurious moment
133  restraintPosition = motion.centreOfRotation();
134 
135  if (motion.report())
136  {
137  Info<< " angle " << theta*sign(a & axis_)
138  << " moment " << restraintMoment
139  << endl;
140  }
141 }
142 
143 
145 (
146  const dictionary& sDoFRBMRDict
147 )
148 {
150 
151  refQ_ = sDoFRBMRCoeffs_.lookupOrDefault<tensor>("referenceOrientation", I);
152 
153  if (mag(mag(refQ_) - sqrt(3.0)) > 1e-9)
154  {
156  << "referenceOrientation " << refQ_ << " is not a rotation tensor. "
157  << "mag(referenceOrientation) - sqrt(3) = "
158  << mag(refQ_) - sqrt(3.0) << nl
159  << exit(FatalError);
160  }
161 
162  axis_ = sDoFRBMRCoeffs_.lookup("axis");
163 
164  scalar magAxis(mag(axis_));
165 
166  if (magAxis > VSMALL)
167  {
168  axis_ /= magAxis;
169  }
170  else
171  {
173  << "axis has zero length"
174  << abort(FatalError);
175  }
176 
177  sDoFRBMRCoeffs_.lookup("stiffness") >> stiffness_;
178  sDoFRBMRCoeffs_.lookup("damping") >> damping_;
179 
180  return true;
181 }
182 
183 
185 (
186  Ostream& os
187 ) const
188 {
189  os.writeKeyword("referenceOrientation")
190  << refQ_ << token::END_STATEMENT << nl;
191 
192  os.writeKeyword("axis")
193  << axis_ << token::END_STATEMENT << nl;
194 
195  os.writeKeyword("stiffness")
196  << stiffness_ << token::END_STATEMENT << nl;
197 
198  os.writeKeyword("damping")
199  << damping_ << token::END_STATEMENT << nl;
200 }
201 
202 // ************************************************************************* //
Six degree of freedom motion for a rigid body.
dimensionedScalar sign(const dimensionedScalar &ds)
virtual void restrain(const sixDoFRigidBodyMotion &motion, vector &restraintPosition, vector &restraintForce, vector &restraintMoment) const
Calculate the restraint position, force and moment.
dimensionedScalar acos(const dimensionedScalar &ds)
bool report() const
Return the report Switch.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
const double e
Elementary charge.
Definition: doubleFloat.H:78
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
const point & centreOfRotation() const
Return the current centre of rotation.
sixDoFRigidBodyMotionRestraints model. Linear axial angular spring.
Base class for defining restraints for sixDoF motions.
virtual bool read(const dictionary &sDoFRBMRCoeff)
Update properties from given dictionary.
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
tensor rotationTensor(const vector &n1, const vector &n2)
Rotational transformation tensor from vector n1 to n2.
Definition: transform.H:47
Macros for easy insertion into run-time selection tables.
bool read(const char *, int32_t &)
Definition: int32IO.C:85
static const Identity< scalar > I
Definition: Identity.H:93
vector omega() const
Return the angular velocity in the global frame.
A class for handling words, derived from string.
Definition: word.H:59
3D tensor transformation operations.
virtual bool read(const dictionary &sDoFRBMRDict)
Update properties from given dictionary.
defineTypeNameAndDebug(linearAxialAngularSpring, 0)
const tensor & orientation() const
Return the orientation tensor, Q.
static const zero Zero
Definition: zero.H:91
errorManip< error > abort(error &err)
Definition: errorManip.H:131
An Ostream is an abstract base class for all output systems (streams, files, token lists...
Definition: Ostream.H:53
static const char nl
Definition: Ostream.H:262
Ostream & writeKeyword(const keyType &)
Write the keyword followed by an appropriate indentation.
Definition: Ostream.C:54
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
addToRunTimeSelectionTable(sixDoFRigidBodyMotionRestraint, linearAxialAngularSpring, dictionary)
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
linearAxialAngularSpring(const word &name, const dictionary &sDoFRBMRDict)
Construct from components.
Namespace for OpenFOAM.