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