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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2017 OpenFOAM Foundation
9  Copyright (C) 2016-2020 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "sixDoFRigidBodyMotion.H"
30 #include "sixDoFSolver.H"
31 #include "septernion.H"
32 
33 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
34 
35 void Foam::sixDoFRigidBodyMotion::applyRestraints()
36 {
37  if (restraints_.empty())
38  {
39  return;
40  }
41 
42  if (Pstream::master())
43  {
44  forAll(restraints_, rI)
45  {
46  if (report_)
47  {
48  Info<< "Restraint " << restraints_[rI].name() << ": ";
49  }
50 
51  // Restraint position
52  point rP = Zero;
53 
54  // Restraint force
55  vector rF = Zero;
56 
57  // Restraint moment
58  vector rM = Zero;
59 
60  // Accumulate the restraints
61  restraints_[rI].restrain(*this, rP, rF, rM);
62 
63  // Update the acceleration
64  a() += rF/mass_;
65 
66  // Moments are returned in global axes, transforming to
67  // body local to add to torque.
68  tau() += Q().T() & (rM + ((rP - centreOfRotation()) ^ rF));
69  }
70  }
71 }
72 
73 
74 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
75 
77 :
78  time_(time),
79  motionState_(),
80  motionState0_(),
81  restraints_(),
82  constraints_(),
83  tConstraints_(tensor::I),
84  rConstraints_(tensor::I),
85  initialCentreOfMass_(Zero),
86  initialCentreOfRotation_(Zero),
87  initialQ_(I),
88  mass_(VSMALL),
89  momentOfInertia_(diagTensor::one*VSMALL),
90  aRelax_(1.0),
91  aDamp_(1.0),
92  report_(false),
93  solver_(nullptr)
94 {}
95 
96 
98 (
99  const dictionary& dict,
100  const dictionary& stateDict,
101  const Time& time
102 )
103 :
104  time_(time),
105  motionState_(stateDict),
106  motionState0_(),
107  restraints_(),
108  constraints_(),
109  tConstraints_(tensor::I),
110  rConstraints_(tensor::I),
111  initialCentreOfMass_
112  (
114  (
115  "initialCentreOfMass",
116  dict.get<vector>("centreOfMass")
117  )
118  ),
119  initialCentreOfRotation_(initialCentreOfMass_),
120  initialQ_
121  (
123  (
124  "initialOrientation",
125  dict.getOrDefault("orientation", tensor::I)
126  )
127  ),
128  mass_(dict.get<scalar>("mass")),
129  momentOfInertia_(dict.get<diagTensor>("momentOfInertia")),
130  aRelax_(dict.getOrDefault<scalar>("accelerationRelaxation", 1)),
131  aDamp_(dict.getOrDefault<scalar>("accelerationDamping", 1)),
132  report_(dict.getOrDefault("report", false)),
133  solver_(sixDoFSolver::New(dict.subDict("solver"), *this))
134 {
135  addRestraints(dict);
136 
137  // Set constraints and initial centre of rotation
138  // if different to the centre of mass
139  addConstraints(dict);
140 
141  // If the centres of mass and rotation are different ...
142  vector R(initialCentreOfMass_ - initialCentreOfRotation_);
143  if (magSqr(R) > VSMALL)
144  {
145  // ... correct the moment of inertia tensor using parallel axes theorem
146  momentOfInertia_ += mass_*diag(I*magSqr(R) - sqr(R));
147 
148  // ... and if the centre of rotation is not specified for motion state
149  // update it
150  if (!stateDict.found("centreOfRotation"))
151  {
152  motionState_.centreOfRotation() = initialCentreOfRotation_;
153  }
154  }
155 
156  // Save the old-time motion state
157  motionState0_ = motionState_;
158 }
159 
160 
162 (
163  const sixDoFRigidBodyMotion& sDoFRBM
164 )
165 :
166  time_(sDoFRBM.time_),
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  solver_(sDoFRBM.solver_.clone())
182 {}
183 
184 
185 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
186 
188 {} // Define here (incomplete type in header)
189 
190 
191 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
192 
194 (
195  const dictionary& dict
196 )
197 {
198  if (dict.found("restraints"))
199  {
200  const dictionary& restraintDict = dict.subDict("restraints");
201 
202  label i = 0;
203 
204  restraints_.setSize(restraintDict.size());
205 
206  for (const entry& dEntry : restraintDict)
207  {
208  if (dEntry.isDict())
209  {
210  restraints_.set
211  (
212  i++,
214  (
215  dEntry.keyword(),
216  dEntry.dict()
217  )
218  );
219  }
220  }
221 
222  restraints_.setSize(i);
223  }
224 }
225 
226 
228 (
229  const dictionary& dict
230 )
231 {
232  if (dict.found("constraints"))
233  {
234  const dictionary& constraintDict = dict.subDict("constraints");
235 
236  label i = 0;
237 
238  constraints_.setSize(constraintDict.size());
239 
240  pointConstraint pct;
241  pointConstraint pcr;
242 
243  for (const entry& dEntry : constraintDict)
244  {
245  if (dEntry.isDict())
246  {
247  constraints_.set
248  (
249  i,
251  (
252  dEntry.keyword(),
253  dEntry.dict(),
254  *this
255  )
256  );
257 
258  constraints_[i].setCentreOfRotation(initialCentreOfRotation_);
259  constraints_[i].constrainTranslation(pct);
260  constraints_[i].constrainRotation(pcr);
261 
262  i++;
263  }
264  }
265 
266  constraints_.setSize(i);
267 
268  tConstraints_ = pct.constraintTransformation();
269  rConstraints_ = pcr.constraintTransformation();
270 
271  Info<< "Translational constraint tensor " << tConstraints_ << nl
272  << "Rotational constraint tensor " << rConstraints_ << endl;
273  }
274 }
275 
276 
277 void Foam::sixDoFRigidBodyMotion::updateAcceleration
278 (
279  const vector& fGlobal,
280  const vector& tauGlobal
281 )
282 {
283  static bool first = true;
284 
285  // Save the previous iteration accelerations for relaxation
286  vector aPrevIter = a();
287  vector tauPrevIter = tau();
288 
289  // Calculate new accelerations
290  a() = fGlobal/mass_;
291  tau() = (Q().T() & tauGlobal);
292  applyRestraints();
293 
294  // Relax accelerations on all but first iteration
295  if (!first)
296  {
297  a() = aRelax_*a() + (1 - aRelax_)*aPrevIter;
298  tau() = aRelax_*tau() + (1 - aRelax_)*tauPrevIter;
299  }
300  else
301  {
302  first = false;
303  }
304 }
305 
306 
308 (
309  bool firstIter,
310  const vector& fGlobal,
311  const vector& tauGlobal,
312  scalar deltaT,
313  scalar deltaT0
314 )
315 {
316  if (Pstream::master())
317  {
318  solver_->solve(firstIter, fGlobal, tauGlobal, deltaT, deltaT0);
319 
320  if (report_)
321  {
322  status();
323  }
324  }
325 
326  Pstream::scatter(motionState_);
327 }
328 
329 
331 {
332  Info<< "6-DoF rigid body motion" << nl
333  << " Centre of rotation: " << centreOfRotation() << nl
334  << " Centre of mass: " << centreOfMass() << nl
335  << " Orientation: " << orientation() << nl
336  << " Linear velocity: " << v() << nl
337  << " Angular velocity: " << omega()
338  << endl;
339 }
340 
341 
343 (
344  const pointField& initialPoints
345 ) const
346 {
347  return
348  (
349  centreOfRotation()
350  + (Q() & initialQ_.T() & (initialPoints - initialCentreOfRotation_))
351  );
352 }
353 
354 
356 (
357  const pointField& initialPoints,
358  const scalarField& scale
359 ) const
360 {
361  // Calculate the transformation septerion from the initial state
362  septernion s
363  (
364  centreOfRotation() - initialCentreOfRotation(),
365  quaternion(Q().T() & initialQ())
366  );
367 
368  tmp<pointField> tpoints(new pointField(initialPoints));
369  pointField& points = tpoints.ref();
370 
371  forAll(points, pointi)
372  {
373  // Move non-stationary points
374  if (scale[pointi] > SMALL)
375  {
376  // Use solid-body motion where scale = 1
377  if (scale[pointi] > 1 - SMALL)
378  {
379  points[pointi] = transform(initialPoints[pointi]);
380  }
381  // Slerp septernion interpolation
382  else
383  {
384  septernion ss(slerp(septernion::I, s, scale[pointi]));
385 
386  points[pointi] =
387  initialCentreOfRotation()
388  + ss.invTransformPoint
389  (
390  initialPoints[pointi]
391  - initialCentreOfRotation()
392  );
393  }
394  }
395  }
396 
397  return tpoints;
398 }
399 
400 
401 // ************************************************************************* //
Foam::entry
A keyword and a list of tokens is an 'entry'.
Definition: entry.H:67
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Foam::Tensor< scalar >
Foam::DiagTensor< scalar >
Foam::Time
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:73
Foam::Tensor::I
static const Tensor I
Definition: Tensor.H:85
Foam::septernion
Septernion class used to perform translations and rotations in 3D space.
Definition: septernion.H:66
s
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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))
Definition: gmvOutputSpray.H:25
Foam::tmp
A class for managing temporary objects.
Definition: PtrList.H:61
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::dictionary::found
bool found(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Search for an entry (const access) with the given keyword.
Definition: dictionaryI.H:87
Foam::slerp
quaternion slerp(const quaternion &qa, const quaternion &qb, const scalar t)
Spherical linear interpolation of quaternions.
Definition: quaternion.C:78
Foam::diag
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
Definition: pointPatchFieldFunctions.H:287
Foam::sixDoFSolver::New
static autoPtr< sixDoFSolver > New(const dictionary &dict, sixDoFRigidBodyMotion &body)
Definition: sixDoFSolverNew.C:34
Foam::pointConstraint
Accumulates point constraints through successive applications of the applyConstraint function.
Definition: pointConstraint.H:60
Foam::one
A class representing the concept of 1 (one) that can be used to avoid manipulating objects known to b...
Definition: one.H:61
Foam::UPstream::master
static bool master(const label communicator=worldComm)
Am I the master process.
Definition: UPstream.H:457
Foam::Pstream::scatter
static void scatter(const List< commsStruct > &comms, T &Value, const int tag, const label comm)
Scatter data. Distribute without modification. Reverse of gather.
Definition: gatherScatter.C:150
Foam::sixDoFRigidBodyMotion::~sixDoFRigidBodyMotion
~sixDoFRigidBodyMotion()
Destructor.
Definition: sixDoFRigidBodyMotion.C:187
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::dictionary::get
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:107
Foam::dictionary::set
entry * set(entry *entryPtr)
Assign a new entry, overwriting any existing entry.
Definition: dictionary.C:780
Foam::transform
dimensionSet transform(const dimensionSet &ds)
Return the argument; transformations do not change the dimensions.
Definition: dimensionSet.C:521
Foam::sixDoFRigidBodyMotion::addRestraints
void addRestraints(const dictionary &dict)
Add restraints to the motion, public to allow external.
Definition: sixDoFRigidBodyMotion.C:194
Foam::sixDoFRigidBodyMotion::sixDoFRigidBodyMotion
sixDoFRigidBodyMotion(const Time &)
Construct null.
Definition: sixDoFRigidBodyMotion.C:76
sixDoFRigidBodyMotion.H
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::quaternion
Quaternion class used to perform rotations in 3D space.
Definition: quaternion.H:56
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
Foam::sixDoFRigidBodyMotionConstraint::New
static autoPtr< sixDoFRigidBodyMotionConstraint > New(const word &name, const dictionary &sDoFRBMCDict, const sixDoFRigidBodyMotion &motion)
Select constructed from the sDoFRBMCDict dictionary and Time.
Definition: sixDoFRigidBodyMotionConstraintNew.C:35
Foam::sixDoFRigidBodyMotion::update
void update(bool firstIter, const vector &fGlobal, const vector &tauGlobal, scalar deltaT, scalar deltaT0)
Symplectic integration of velocities, orientation and position.
Definition: sixDoFRigidBodyMotion.C:308
R
#define R(A, B, C, D, E, F, K, M)
Foam::sixDoFRigidBodyMotion::centreOfRotation
const point & centreOfRotation() const
Return the current centre of rotation.
Definition: sixDoFRigidBodyMotionI.H:231
Foam::tmp::ref
T & ref() const
Definition: tmpI.H:227
septernion.H
Foam::Field< vector >
Foam::sixDoFRigidBodyMotion::addConstraints
void addConstraints(const dictionary &dict)
Add restraints to the motion, public to allow external.
Definition: sixDoFRigidBodyMotion.C:228
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::dictionary::subDict
const dictionary & subDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary.
Definition: dictionary.C:460
Foam::T
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
Definition: FieldFieldFunctions.C:58
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:51
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::sixDoFRigidBodyMotion
Six degree of freedom motion for a rigid body.
Definition: sixDoFRigidBodyMotion.H:69
Foam::sixDoFRigidBodyMotion::transform
point transform(const point &initialPoints) const
Transform the given initial state point by the current motion.
Definition: sixDoFRigidBodyMotionI.H:301
Foam::nl
constexpr char nl
Definition: Ostream.H:404
Foam::Vector< scalar >
points
const pointField & points
Definition: gmvOutputHeader.H:1
sixDoFSolver.H
Foam::Tensor::T
Tensor< Cmpt > T() const
Return non-Hermitian transpose.
Definition: TensorI.H:504
Foam::septernion::I
static const septernion I
Definition: septernion.H:84
Foam::pointConstraint::constraintTransformation
tensor constraintTransformation() const
Return the accumulated constraint transformation tensor.
Definition: pointConstraintI.H:120
Foam::septernion::invTransformPoint
vector invTransformPoint(const vector &v) const
Inverse Transform the given coordinate point.
Definition: septernionI.H:98
Foam::point
vector point
Point is a vector.
Definition: point.H:43
Foam::dictionary::getOrDefault
T getOrDefault(const word &keyword, const T &deflt, enum keyType::option matchOpt=keyType::REGEX) const
Definition: dictionaryTemplates.C:148
Foam::sixDoFRigidBodyMotion::status
void status() const
Report the status of the motion.
Definition: sixDoFRigidBodyMotion.C:330
Foam::sixDoFRigidBodyMotionRestraint::New
static autoPtr< sixDoFRigidBodyMotionRestraint > New(const word &name, const dictionary &sDoFRBMRDict)
Select constructed from the sDoFRBMRDict dictionary and Time.
Definition: sixDoFRigidBodyMotionRestraintNew.C:35
Foam::I
static const Identity< scalar > I
Definition: Identity.H:95