rigidBodyMotion.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) 2016-2017 OpenFOAM Foundation
9  Copyright (C) 2020-2021 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 "rigidBodyMotion.H"
30 #include "rigidBodySolver.H"
31 #include "septernion.H"
32 
33 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
34 
35 void Foam::RBD::rigidBodyMotion::initialize()
36 {
37  // Calculate the initial body-state
38  forwardDynamicsCorrection(rigidBodyModelState(*this));
39  X00_ = X0_;
40 
41  // Update the body-state to correspond to the current joint-state
42  forwardDynamicsCorrection(motionState_);
43 }
44 
45 
46 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
47 
48 Foam::RBD::rigidBodyMotion::rigidBodyMotion(const Time& time)
49 :
50  rigidBodyModel(time),
51  motionState_(*this),
52  motionState0_(*this),
53  aRelax_(nullptr),
54  aDamp_(1.0),
55  report_(false),
56  solver_(nullptr)
57 {}
58 
59 Foam::RBD::rigidBodyMotion::rigidBodyMotion
60 (
61  const Time& time,
62  const dictionary& dict
63 )
64 :
65  rigidBodyModel(time, dict),
66  motionState_(*this, dict),
67  motionState0_(motionState_),
68  X00_(X0_.size()),
69  aRelax_
70  (
72  (
73  "accelerationRelaxation",
74  dict,
75  word::null,
76  &time
77  )
78  ),
79  aDamp_(dict.getOrDefault<scalar>("accelerationDamping", 1)),
80  report_(dict.getOrDefault<Switch>("report", false)),
81  solver_(rigidBodySolver::New(*this, dict.subDict("solver")))
82 {
83  if (dict.found("g"))
84  {
85  g() = dict.get<vector>("g");
86  }
87 
88  initialize();
89 }
90 
91 
92 Foam::RBD::rigidBodyMotion::rigidBodyMotion
93 (
94  const Time& time,
95  const dictionary& dict,
96  const dictionary& stateDict
97 )
98 :
99  rigidBodyModel(time, dict),
100  motionState_(*this, stateDict),
101  motionState0_(motionState_),
102  X00_(X0_.size()),
103  aRelax_
104  (
106  (
107  "accelerationRelaxation",
108  dict,
109  word::null,
110  &time
111  )
112  ),
113  aDamp_(dict.getOrDefault<scalar>("accelerationDamping", 1)),
114  report_(dict.getOrDefault<Switch>("report", false)),
115  solver_(rigidBodySolver::New(*this, dict.subDict("solver")))
116 {
117  if (dict.found("g"))
118  {
119  g() = dict.get<vector>("g");
120  }
121 
122  initialize();
123 }
124 
125 
126 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
127 
129 {}
130 
131 
132 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
133 
135 (
136  const label bodyId
137 ) const
138 {
139  if (merged(bodyId))
140  {
141  const subBody& mBody = mergedBody(bodyId);
142  return mBody.masterXT() & X00_[mBody.masterID()];
143  }
144  else
145  {
146  return X00_[bodyId];
147  }
148 }
149 
150 
152 (
153  rigidBodyModelState& state,
154  const scalarField& tau,
155  const Field<spatialVector>& fx
156 ) const
157 {
158  scalarField qDdotPrev = state.qDdot();
159  rigidBodyModel::forwardDynamics(state, tau, fx);
160 
161  scalar aRelax = 1;
162  if (aRelax_)
163  {
164  aRelax = aRelax_->value(motionState_.t());
165  }
166 
167  state.qDdot() = aDamp_*(aRelax*state.qDdot() + (1 - aRelax)*qDdotPrev);
168 }
169 
170 
172 (
173  const scalar t,
174  const scalar deltaT,
175  const scalarField& tau,
176  const Field<spatialVector>& fx
177 )
178 {
179  motionState_.t() = t;
180  motionState_.deltaT() = deltaT;
181 
182  if (motionState0_.deltaT() < SMALL)
183  {
184  motionState0_.t() = t;
185  motionState0_.deltaT() = deltaT;
186  }
187 
188  if (Pstream::master())
189  {
190  solver_->solve(tau, fx);
191  }
192 
193  Pstream::scatter(motionState_);
194 
195  // Update the body-state to correspond to the current joint-state
196  forwardDynamicsCorrection(motionState_);
197 }
198 
199 
200 void Foam::RBD::rigidBodyMotion::status(const label bodyID) const
201 {
202  const spatialTransform CofR(X0(bodyID));
203  const spatialVector vCofR(v(bodyID, Zero));
204 
205  Info<< "Rigid-body motion of the " << name(bodyID) << nl
206  << " Centre of rotation: " << CofR.r() << nl
207  << " Orientation: " << CofR.E() << nl
208  << " Linear velocity: " << vCofR.l() << nl
209  << " Angular velocity: " << vCofR.w()
210  << endl;
211 }
212 
213 
214 const Foam::vector Foam::RBD::rigidBodyMotion::vCofR(const label bodyID) const
215 {
216  const spatialVector velCofR(v(bodyID, Zero));
217 
218  return velCofR.l();
219 }
220 
221 
222 const Foam::vector Foam::RBD::rigidBodyMotion::cCofR(const label bodyID) const
223 {
224  const spatialTransform CofR(X0(bodyID));
225 
226  return CofR.r();
227 }
228 
229 
230 
232 (
233  const label bodyID,
234  const pointField& initialPoints
235 ) const
236 {
237  // Calculate the transform from the initial state in the global frame
238  // to the current state in the global frame
239  spatialTransform X(X0(bodyID).inv() & X00(bodyID));
240 
241  tmp<pointField> tpoints(new pointField(initialPoints.size()));
242  pointField& points = tpoints.ref();
243 
244  forAll(points, i)
245  {
246  points[i] = X.transformPoint(initialPoints[i]);
247  }
248 
249  return tpoints;
250 }
251 
252 
254 (
255  const label bodyID,
256  const scalarField& weight,
257  const pointField& initialPoints
258 ) const
259 {
260  // Calculate the transform from the initial state in the global frame
261  // to the current state in the global frame
262  spatialTransform X(X0(bodyID).inv() & X00(bodyID));
263 
264  // Calculate the septernion equivalent of the transformation for 'slerp'
265  // interpolation
266  septernion s(X);
267 
268  tmp<pointField> tpoints(new pointField(initialPoints));
269  pointField& points = tpoints.ref();
270 
271  forAll(points, i)
272  {
273  // Move non-stationary points
274  if (weight[i] > SMALL)
275  {
276  // Use solid-body motion where weight = 1
277  if (weight[i] > 1 - SMALL)
278  {
279  points[i] = X.transformPoint(initialPoints[i]);
280  }
281  // Slerp septernion interpolation
282  else
283  {
284  points[i] =
285  slerp(septernion::I, s, weight[i])
286  .transformPoint(initialPoints[i]);
287  }
288  }
289  }
290 
291  return tpoints;
292 }
293 
294 
296 (
297  const labelList& bodyIDs,
298  const List<const scalarField*>& weights,
299  const pointField& initialPoints
300 ) const
301 {
302  List<septernion> ss(bodyIDs.size() + 1);
303  ss[bodyIDs.size()] = septernion::I;
304 
305  forAll(bodyIDs, bi)
306  {
307  const label bodyID = bodyIDs[bi];
308 
309  // Calculate the transform from the initial state in the global frame
310  // to the current state in the global frame
311  spatialTransform X(X0(bodyID).inv() & X00(bodyID));
312 
313  // Calculate the septernion equivalent of the transformation
314  ss[bi] = septernion(X);
315  }
316 
317  tmp<pointField> tpoints(new pointField(initialPoints));
318  pointField& points = tpoints.ref();
319 
320  List<scalar> w(ss.size());
321 
322  forAll(points, i)
323  {
324  // Initialize to 1 for the far-field weight
325  scalar sum1mw = 1;
326 
327  forAll(bodyIDs, bi)
328  {
329  w[bi] = (*(weights[bi]))[i];
330  sum1mw += w[bi]/(1 + SMALL - w[bi]);
331  }
332 
333  // Calculate the limiter for wi/(1 - wi) to ensure the sum(wi) = 1
334  scalar lambda = 1/sum1mw;
335 
336  // Limit wi/(1 - wi) and sum the resulting wi
337  scalar sumw = 0;
338  forAll(bodyIDs, bi)
339  {
340  w[bi] = lambda*w[bi]/(1 + SMALL - w[bi]);
341  sumw += w[bi];
342  }
343 
344  // Calculate the weight for the stationary far-field
345  w[bodyIDs.size()] = 1 - sumw;
346 
347  points[i] = average(ss, w).transformPoint(initialPoints[i]);
348  }
349 
350  return tpoints;
351 }
352 
353 
354 // ************************************************************************* //
Foam::RBD::rigidBodyMotion::solve
void solve(const scalar t, const scalar deltaT, const scalarField &tau, const Field< spatialVector > &fx)
Integrate velocities, orientation and position.
Definition: rigidBodyMotion.C:172
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Foam::RBD::rigidBodyMotion::cCofR
const vector cCofR(const label bodyID) const
Report CofR of the given body.
Definition: rigidBodyMotion.C:222
rigidBodySolver.H
Foam::Switch
A simple wrapper around bool so that it can be read as a word: true/false, on/off,...
Definition: Switch.H:77
Foam::Time
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:73
Foam::RBD::rigidBodyModelState
Holds the motion state of rigid-body model.
Definition: rigidBodyModelState.H:67
Foam::RBD::rigidBodyModelState::qDdot
const scalarField & qDdot() const
Return access to the joint acceleration.
Definition: rigidBodyModelStateI.H:42
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::spatialTransform
Compact representation of the Plücker spatial transformation tensor in terms of the rotation tensor E...
Definition: spatialTransform.H:70
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
X0
scalarList X0(nSpecie, Zero)
Foam::slerp
quaternion slerp(const quaternion &qa, const quaternion &qb, const scalar t)
Spherical linear interpolation of quaternions.
Definition: quaternion.C:78
Foam::RBD::rigidBodyMotion::status
void status(const label bodyID) const
Report the status of the motion of the given body.
Definition: rigidBodyMotion.C:200
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::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
Foam::RBD::subBody::masterID
label masterID() const
Return the master body Id.
Definition: subBodyI.H:71
Foam::Function1
Top level data entry class for use in dictionaries. Provides a mechanism to specify a variable as a c...
Definition: propellerInfo.H:291
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::RBD::rigidBodyModel::X0_
DynamicList< spatialTransform > X0_
Transform for external forces to the bodies reference frame.
Definition: rigidBodyModel.H:154
Foam::RBD::rigidBodySolver::New
static autoPtr< rigidBodySolver > New(rigidBodyMotion &body, const dictionary &dict)
Definition: rigidBodySolverNew.C:34
Foam::RBD::rigidBodyModel::forwardDynamicsCorrection
void forwardDynamicsCorrection(const rigidBodyModelState &state) const
Correct the velocity and acceleration of the bodies in the model.
Definition: forwardDynamics.C:210
Foam::tmp::ref
T & ref() const
Definition: tmpI.H:227
septernion.H
Foam::Field< scalar >
Foam::spatialTransform::r
const vector & r() const
Return the translation vector.
Definition: spatialTransformI.H:89
Foam::inv
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
Definition: dimensionedSphericalTensor.C:73
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::SpatialVector::l
Vector< Cmpt > l() const
Return the linear part of the spatial vector as a vector.
Definition: SpatialVectorI.H:195
Foam::SpatialVector< scalar >
Foam::RBD::subBody
Definition: subBody.H:57
Foam::spatialTransform::E
const tensor & E() const
Return the rotation tensor.
Definition: spatialTransformI.H:79
Foam::RBD::rigidBodyMotion::vCofR
const vector vCofR(const label bodyID) const
Report linear velocity of the given body.
Definition: rigidBodyMotion.C:214
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::RBD::rigidBodyModel
Basic rigid-body model representing a system of rigid-bodies connected by 1-6 DoF joints.
Definition: rigidBodyModel.H:83
rigidBodyMotion.H
Foam::spatialTransform::transformPoint
vector transformPoint(const vector &p) const
Transform position p.
Definition: spatialTransformI.H:173
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
lambda
dimensionedScalar lambda("lambda", dimTime/sqr(dimLength), laminarTransport)
Foam::RBD::rigidBodyMotion::X00
spatialTransform X00(const label bodyId) const
Return the initial transform to the global frame for the.
Definition: rigidBodyMotion.C:135
g
const uniformDimensionedVectorField & g
Definition: createFluidFields.H:26
Foam::RBD::rigidBodyModel::forwardDynamics
void forwardDynamics(rigidBodyModelState &state, const scalarField &tau, const Field< spatialVector > &fx) const
Calculate the joint acceleration qDdot from the joint state q,.
Definition: forwardDynamics.C:57
Foam::nl
constexpr char nl
Definition: Ostream.H:404
Foam::Vector< scalar >
Foam::List< label >
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::word::null
static const word null
An empty word.
Definition: word.H:80
Foam::RBD::rigidBodyMotion::~rigidBodyMotion
~rigidBodyMotion()
Destructor.
Definition: rigidBodyMotion.C:128
Foam::septernion::I
static const septernion I
Definition: septernion.H:84
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
Foam::RBD::subBody::masterXT
const spatialTransform & masterXT() const
Return the transform with respect to the master body.
Definition: subBodyI.H:77
Foam::RBD::rigidBodyMotion::transformPoints
tmp< pointField > transformPoints(const label bodyID, const pointField &initialPoints) const
Transform the given initial pointField of the specified body.
Definition: rigidBodyMotion.C:232
Foam::RBD::rigidBodyMotion::forwardDynamics
void forwardDynamics(rigidBodyModelState &state, const scalarField &tau, const Field< spatialVector > &fx) const
Calculate and optionally relax the joint acceleration qDdot from.
Definition: rigidBodyMotion.C:152
Foam::SpatialVector::w
Vector< Cmpt > w() const
Return the angular part of the spatial vector as a vector.
Definition: SpatialVectorI.H:189
Foam::average
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
Definition: DimensionedFieldFunctions.C:328