FreeFOAM The Cross-Platform CFD Toolkit
tabulatedAxialAngularSpring.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) 1991-2010 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 
29 #include <OpenFOAM/transform.H>
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 namespace sixDoFRigidBodyMotionRestraints
37 {
38  defineTypeNameAndDebug(tabulatedAxialAngularSpring, 0);
40  (
41  sixDoFRigidBodyMotionRestraint,
42  tabulatedAxialAngularSpring,
43  dictionary
44  );
45 };
46 };
47 
48 
49 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
50 
53 (
54  const dictionary& sDoFRBMRDict
55 )
56 :
57  sixDoFRigidBodyMotionRestraint(sDoFRBMRDict),
58  refQ_(),
59  axis_(),
60  moment_(),
61  convertToDegrees_(),
62  damping_()
63 {
64  read(sDoFRBMRDict);
65 }
66 
67 
68 // * * * * * * * * * * * * * * * * Destructors * * * * * * * * * * * * * * * //
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 
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  vector oldDir = refQ_ & refDir;
99 
100  vector newDir = motion.orientation() & refDir;
101  }
102 
103  // Removing any axis component from oldDir and newDir and normalising
104  oldDir -= (axis_ & oldDir)*axis_;
105  oldDir /= (mag(oldDir) + VSMALL);
106 
107  newDir -= (axis_ & newDir)*axis_;
108  newDir /= (mag(newDir) + VSMALL);
109 
110  scalar theta = mag(acos(min(oldDir & newDir, 1.0)));
111 
112  // Determining the sign of theta by comparing the cross product of
113  // the direction vectors with the axis
114  theta *= sign((oldDir ^ newDir) & axis_);
115 
116  scalar moment;
117 
118  if (convertToDegrees_)
119  {
120  moment = moment_(theta*180.0/mathematicalConstant::pi);
121  }
122  else
123  {
124  moment = moment_(theta);
125  }
126 
127  // Damping of along axis angular velocity only
128  restraintMoment = moment*axis_ - damping_*(motion.omega() & axis_)*axis_;
129 
130  restraintForce = vector::zero;
131 
132  // Not needed to be altered as restraintForce is zero, but set to
133  // centreOfMass to be sure of no spurious moment
134  restraintPosition = motion.centreOfMass();
135 
136  if (motion.report())
137  {
138  Info<< " angle " << theta
139  << " force " << restraintForce
140  << " moment " << restraintMoment
141  << endl;
142  }
143 }
144 
145 
147 (
148  const dictionary& sDoFRBMRDict
149 )
150 {
152 
153  refQ_ = sDoFRBMRCoeffs_.lookupOrDefault<tensor>("referenceOrientation", I);
154 
155  if (mag(mag(refQ_) - sqrt(3.0)) > 1e-9)
156  {
158  (
159  "Foam::sixDoFRigidBodyMotionRestraints::"
160  "tabulatedAxialAngularSpring::read"
161  "("
162  "const dictionary& sDoFRBMRDict"
163  ")"
164  )
165  << "referenceOrientation " << refQ_ << " is not a rotation tensor. "
166  << "mag(referenceOrientation) - sqrt(3) = "
167  << mag(refQ_) - sqrt(3.0) << nl
168  << exit(FatalError);
169  }
170 
171  axis_ = sDoFRBMRCoeffs_.lookup("axis");
172 
173  scalar magAxis(mag(axis_));
174 
175  if (magAxis > VSMALL)
176  {
177  axis_ /= magAxis;
178  }
179  else
180  {
182  (
183  "Foam::sixDoFRigidBodyMotionRestraints::"
184  "tabulatedAxialAngularSpring::read"
185  "("
186  "const dictionary& sDoFRBMCDict"
187  ")"
188  )
189  << "axis has zero length"
190  << abort(FatalError);
191  }
192 
193  moment_ = interpolationTable<scalar>(sDoFRBMRCoeffs_);
194 
195  word angleFormat = sDoFRBMRCoeffs_.lookup("angleFormat");
196 
197  if (angleFormat == "degrees" || angleFormat == "degree")
198  {
199  convertToDegrees_ = true;
200  }
201  else if (angleFormat == "radians" || angleFormat == "radian")
202  {
203  convertToDegrees_ = false;
204  }
205  else
206  {
208  (
209  "Foam::sixDoFRigidBodyMotionRestraints::"
210  "tabulatedAxialAngularSpring::read"
211  "("
212  "const dictionary& sDoFRBMCDict"
213  ")"
214  )
215  << "angleFormat must be degree, degrees, radian or radians"
216  << abort(FatalError);
217  }
218 
219  sDoFRBMRCoeffs_.lookup("damping") >> damping_;
220 
221  return true;
222 }
223 
224 
226 (
227  Ostream& os
228 ) const
229 {
230  os.writeKeyword("referenceOrientation")
231  << refQ_ << token::END_STATEMENT << nl;
232 
233  os.writeKeyword("axis")
234  << axis_ << token::END_STATEMENT << nl;
235 
236  moment_.write(os);
237 
238  os.writeKeyword("angleFormat");
239 
240  if (convertToDegrees_)
241  {
242  os << "degrees" << token::END_STATEMENT << nl;
243  }
244  else
245  {
246  os << "radians" << token::END_STATEMENT << nl;
247  }
248 
249  os.writeKeyword("damping")
250  << damping_ << token::END_STATEMENT << nl;
251 }
252 
253 
254 // ************************ vim: set sw=4 sts=4 et: ************************ //