sixDoFRigidBodyMotion.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-2015 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 
26 #include "sixDoFRigidBodyMotion.H"
27 #include "sixDoFSolver.H"
28 #include "septernion.H"
29 
30 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
31 
32 void Foam::sixDoFRigidBodyMotion::applyRestraints()
33 {
34  if (restraints_.empty())
35  {
36  return;
37  }
38 
39  if (Pstream::master())
40  {
41  forAll(restraints_, rI)
42  {
43  if (report_)
44  {
45  Info<< "Restraint " << restraints_[rI].name() << ": ";
46  }
47 
48  // Restraint position
49  point rP = vector::zero;
50 
51  // Restraint force
52  vector rF = vector::zero;
53 
54  // Restraint moment
55  vector rM = vector::zero;
56 
57  // Accumulate the restraints
58  restraints_[rI].restrain(*this, rP, rF, rM);
59 
60  // Update the acceleration
61  a() += rF/mass_;
62 
63  // Moments are returned in global axes, transforming to
64  // body local to add to torque.
65  tau() += Q().T() & (rM + ((rP - centreOfRotation()) ^ rF));
66  }
67  }
68 }
69 
70 
71 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
72 
74 :
75  motionState_(),
76  motionState0_(),
77  restraints_(),
78  constraints_(),
79  tConstraints_(tensor::I),
80  rConstraints_(tensor::I),
81  initialCentreOfMass_(vector::zero),
82  initialCentreOfRotation_(vector::zero),
83  initialQ_(I),
84  mass_(VSMALL),
85  momentOfInertia_(diagTensor::one*VSMALL),
86  aRelax_(1.0),
87  aDamp_(1.0),
88  report_(false),
89  solver_(NULL)
90 {}
91 
92 
94 (
95  const dictionary& dict,
96  const dictionary& stateDict
97 )
98 :
99  motionState_(stateDict),
100  motionState0_(),
101  restraints_(),
102  constraints_(),
103  tConstraints_(tensor::I),
104  rConstraints_(tensor::I),
105  initialCentreOfMass_
106  (
107  dict.lookupOrDefault
108  (
109  "initialCentreOfMass",
110  vector(dict.lookup("centreOfMass"))
111  )
112  ),
113  initialCentreOfRotation_(initialCentreOfMass_),
114  initialQ_
115  (
116  dict.lookupOrDefault
117  (
118  "initialOrientation",
119  dict.lookupOrDefault("orientation", tensor::I)
120  )
121  ),
122  mass_(readScalar(dict.lookup("mass"))),
123  momentOfInertia_(dict.lookup("momentOfInertia")),
124  aRelax_(dict.lookupOrDefault<scalar>("accelerationRelaxation", 1.0)),
125  aDamp_(dict.lookupOrDefault<scalar>("accelerationDamping", 1.0)),
126  report_(dict.lookupOrDefault<Switch>("report", false)),
127  solver_(sixDoFSolver::New(dict.subDict("solver"), *this))
128 {
129  addRestraints(dict);
130 
131  // Set constraints and initial centre of rotation
132  // if different to the centre of mass
133  addConstraints(dict);
134 
135  // If the centres of mass and rotation are different ...
136  vector R(initialCentreOfMass_ - initialCentreOfRotation_);
137  if (magSqr(R) > VSMALL)
138  {
139  // ... correct the moment of inertia tensor using parallel axes theorem
140  momentOfInertia_ += mass_*diag(I*magSqr(R) - sqr(R));
141 
142  // ... and if the centre of rotation is not specified for motion state
143  // update it
144  if (!stateDict.found("centreOfRotation"))
145  {
146  motionState_.centreOfRotation() = initialCentreOfRotation_;
147  }
148  }
149 
150  // Save the old-time motion state
151  motionState0_ = motionState_;
152 }
153 
154 
156 (
157  const sixDoFRigidBodyMotion& sDoFRBM
158 )
159 :
160  motionState_(sDoFRBM.motionState_),
161  motionState0_(sDoFRBM.motionState0_),
162  restraints_(sDoFRBM.restraints_),
163  constraints_(sDoFRBM.constraints_),
164  tConstraints_(sDoFRBM.tConstraints_),
165  rConstraints_(sDoFRBM.rConstraints_),
166  initialCentreOfMass_(sDoFRBM.initialCentreOfMass_),
167  initialCentreOfRotation_(sDoFRBM.initialCentreOfRotation_),
168  initialQ_(sDoFRBM.initialQ_),
169  mass_(sDoFRBM.mass_),
170  momentOfInertia_(sDoFRBM.momentOfInertia_),
171  aRelax_(sDoFRBM.aRelax_),
172  aDamp_(sDoFRBM.aDamp_),
173  report_(sDoFRBM.report_)
174 {}
175 
176 
177 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
178 
180 {}
181 
182 
183 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
184 
186 (
187  const dictionary& dict
188 )
189 {
190  if (dict.found("restraints"))
191  {
192  const dictionary& restraintDict = dict.subDict("restraints");
193 
194  label i = 0;
195 
196  restraints_.setSize(restraintDict.size());
197 
198  forAllConstIter(IDLList<entry>, restraintDict, iter)
199  {
200  if (iter().isDict())
201  {
202  restraints_.set
203  (
204  i++,
206  (
207  iter().keyword(),
208  iter().dict()
209  )
210  );
211  }
212  }
213 
214  restraints_.setSize(i);
215  }
216 }
217 
218 
220 (
221  const dictionary& dict
222 )
223 {
224  if (dict.found("constraints"))
225  {
226  const dictionary& constraintDict = dict.subDict("constraints");
227 
228  label i = 0;
229 
230  constraints_.setSize(constraintDict.size());
231 
232  pointConstraint pct;
233  pointConstraint pcr;
234 
235  forAllConstIter(IDLList<entry>, constraintDict, iter)
236  {
237  if (iter().isDict())
238  {
239  constraints_.set
240  (
241  i,
243  (
244  iter().keyword(),
245  iter().dict(),
246  *this
247  )
248  );
249 
250  constraints_[i].setCentreOfRotation(initialCentreOfRotation_);
251  constraints_[i].constrainTranslation(pct);
252  constraints_[i].constrainRotation(pcr);
253 
254  i++;
255  }
256  }
257 
258  constraints_.setSize(i);
259 
260  tConstraints_ = pct.constraintTransformation();
261  rConstraints_ = pcr.constraintTransformation();
262 
263  Info<< "Translational constraint tensor " << tConstraints_ << nl
264  << "Rotational constraint tensor " << rConstraints_ << endl;
265  }
266 }
267 
268 
269 void Foam::sixDoFRigidBodyMotion::updateAcceleration
270 (
271  const vector& fGlobal,
272  const vector& tauGlobal
273 )
274 {
275  static bool first = false;
276 
277  // Save the previous iteration accelerations for relaxation
278  vector aPrevIter = a();
279  vector tauPrevIter = tau();
280 
281  // Calculate new accelerations
282  a() = fGlobal/mass_;
283  tau() = (Q().T() & tauGlobal);
284  applyRestraints();
285 
286  // Relax accelerations on all but first iteration
287  if (!first)
288  {
289  a() = aRelax_*a() + (1 - aRelax_)*aPrevIter;
290  tau() = aRelax_*tau() + (1 - aRelax_)*tauPrevIter;
291  }
292 
293  first = false;
294 }
295 
296 
298 (
299  bool firstIter,
300  const vector& fGlobal,
301  const vector& tauGlobal,
302  scalar deltaT,
303  scalar deltaT0
304 )
305 {
306  if (Pstream::master())
307  {
308  solver_->solve(firstIter, fGlobal, tauGlobal, deltaT, deltaT0);
309 
310  if (report_)
311  {
312  status();
313  }
314  }
315 
316  Pstream::scatter(motionState_);
317 }
318 
319 
321 {
322  Info<< "6-DoF rigid body motion" << nl
323  << " Centre of rotation: " << centreOfRotation() << nl
324  << " Centre of mass: " << centreOfMass() << nl
325  << " Orientation: " << orientation() << nl
326  << " Linear velocity: " << v() << nl
327  << " Angular velocity: " << omega()
328  << endl;
329 }
330 
331 
333 (
334  const pointField& initialPoints
335 ) const
336 {
337  return
338  (
340  + (Q() & initialQ_.T() & (initialPoints - initialCentreOfRotation_))
341  );
342 }
343 
344 
346 (
347  const pointField& initialPoints,
348  const scalarField& scale
349 ) const
350 {
351  // Calculate the transformation septerion from the initial state
352  septernion s
353  (
354  centreOfRotation() - initialCentreOfRotation(),
355  quaternion(Q() & initialQ().T())
356  );
357 
358  tmp<pointField> tpoints(new pointField(initialPoints));
359  pointField& points = tpoints();
360 
361  forAll(points, pointi)
362  {
363  // Move non-stationary points
364  if (scale[pointi] > SMALL)
365  {
366  // Use solid-body motion where scale = 1
367  if (scale[pointi] > 1 - SMALL)
368  {
369  points[pointi] = transform(initialPoints[pointi]);
370  }
371  // Slerp septernion interpolation
372  else
373  {
374  septernion ss(slerp(septernion::I, s, scale[pointi]));
375 
376  points[pointi] =
377  initialCentreOfRotation()
378  + ss.transform
379  (
380  initialPoints[pointi]
381  - initialCentreOfRotation()
382  );
383  }
384  }
385  }
386 
387  return tpoints;
388 }
389 
390 
391 // ************************************************************************* //
bool found(const word &, bool recursive=false, bool patternMatch=true) const
Search dictionary for given keyword.
Definition: dictionary.C:306
const pointField & points
vector transform(const vector &v) const
Transform the given vector.
Definition: septernionI.H:76
static autoPtr< sixDoFSolver > New(const dictionary &dict, sixDoFRigidBodyMotion &body)
A class representing the concept of 1 (scalar(1.0)) used to avoid unnecessary manipulations for objec...
Definition: one.H:45
vector point
Point is a vector.
Definition: point.H:41
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject( name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE ))
void status() const
Report the status of the motion.
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
dimensioned< scalar > magSqr(const dimensioned< Type > &)
#define R(A, B, C, D, E, F, K, M)
const tensor & orientation() const
Return the orientation tensor, Q.
static const sphericalTensor I(1)
Quaternion class used to perform rotations in 3D space.
Definition: quaternion.H:60
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
A simple wrapper around bool so that it can be read as a word: true/false, on/off, yes/no, y/n, t/f, or none.
Definition: Switch.H:60
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
Accumulates point constraints through successive applications of the applyConstraint function...
messageStream Info
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:638
label size() const
Return number of elements in list.
Definition: DLListBaseI.H:77
void addRestraints(const dictionary &dict)
Add restraints to the motion, public to allow external.
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
point transform(const point &initialPoints) const
Transform the given initial state point by the current motion.
Septernion class used to perform translations and rotations in 3D space.
Definition: septernion.H:64
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:39
dictionary dict
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
static autoPtr< sixDoFRigidBodyMotionConstraint > New(const word &name, const dictionary &sDoFRBMCDict, const sixDoFRigidBodyMotion &motion)
Select constructed from the sDoFRBMCDict dictionary and Time.
void update(bool firstIter, const vector &fGlobal, const vector &tauGlobal, scalar deltaT, scalar deltaT0)
Symplectic integration of velocities, orientation and position.
#define forAll(list, i)
Definition: UList.H:421
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
A class representing the concept of 0 used to avoid unnecessary manipulations for objects that are kn...
Definition: zero.H:47
const point & centreOfRotation() const
Return access to the centre of mass.
static void scatter(const List< commsStruct > &comms, T &Value, const int tag, const label comm)
Scatter data. Distribute without modification. Reverse of gather.
void addConstraints(const dictionary &dict)
Add restraints to the motion, public to allow external.
tensor constraintTransformation() const
Return the accumulated constraint transformation tensor.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:452
static autoPtr< sixDoFRigidBodyMotionRestraint > New(const word &name, const dictionary &sDoFRBMRDict)
Select constructed from the sDoFRBMRDict dictionary and Time.
T lookupOrDefault(const word &, const T &, bool recursive=false, bool patternMatch=true) const
Find and return a T,.
vector omega() const
Return the angular velocity in the global frame.
static const Vector zero
Definition: Vector.H:80
Tensor< Cmpt > T() const
Transpose.
Definition: TensorI.H:286
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:398
static const septernion I
Definition: septernion.H:82
Six degree of freedom motion for a rigid body.
sixDoFRigidBodyMotion()
Construct null.
static const Tensor I
Definition: Tensor.H:84
dimensionedSymmTensor sqr(const dimensionedVector &dv)
quaternion slerp(const quaternion &qa, const quaternion &qb, const scalar t)
Spherical linear interpolation of quaternions.
Definition: quaternion.C:55
This function object calculates and outputs the second invariant of the velocity gradient tensor [1/s...
Definition: Q.H:66
const point & centreOfRotation() const
Return the current centre of rotation.
point centreOfMass() const
Return the current centre of mass.
A class for managing temporary objects.
Definition: PtrList.H:118