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 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_(1.0),
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_(dict.getOrDefault<scalar>("accelerationRelaxation", 1)),
70  aDamp_(dict.getOrDefault<scalar>("accelerationDamping", 1)),
71  report_(dict.getOrDefault<Switch>("report", false)),
72  solver_(rigidBodySolver::New(*this, dict.subDict("solver")))
73 {
74  if (dict.found("g"))
75  {
76  g() = dict.get<vector>("g");
77  }
78 
79  initialize();
80 }
81 
82 
83 Foam::RBD::rigidBodyMotion::rigidBodyMotion
84 (
85  const Time& time,
86  const dictionary& dict,
87  const dictionary& stateDict
88 )
89 :
90  rigidBodyModel(time, dict),
91  motionState_(*this, stateDict),
92  motionState0_(motionState_),
93  X00_(X0_.size()),
94  aRelax_(dict.getOrDefault<scalar>("accelerationRelaxation", 1)),
95  aDamp_(dict.getOrDefault<scalar>("accelerationDamping", 1)),
96  report_(dict.getOrDefault<Switch>("report", false)),
97  solver_(rigidBodySolver::New(*this, dict.subDict("solver")))
98 {
99  if (dict.found("g"))
100  {
101  g() = dict.get<vector>("g");
102  }
103 
104  initialize();
105 }
106 
107 
108 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
109 
111 {}
112 
113 
114 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
115 
117 (
118  const label bodyId
119 ) const
120 {
121  if (merged(bodyId))
122  {
123  const subBody& mBody = mergedBody(bodyId);
124  return mBody.masterXT() & X00_[mBody.masterID()];
125  }
126  else
127  {
128  return X00_[bodyId];
129  }
130 }
131 
132 
134 (
135  rigidBodyModelState& state,
136  const scalarField& tau,
137  const Field<spatialVector>& fx
138 ) const
139 {
140  scalarField qDdotPrev = state.qDdot();
141  rigidBodyModel::forwardDynamics(state, tau, fx);
142  state.qDdot() = aDamp_*(aRelax_*state.qDdot() + (1 - aRelax_)*qDdotPrev);
143 }
144 
145 
147 (
148  const scalar t,
149  const scalar deltaT,
150  const scalarField& tau,
151  const Field<spatialVector>& fx
152 )
153 {
154  motionState_.t() = t;
155  motionState_.deltaT() = deltaT;
156 
157  if (motionState0_.deltaT() < SMALL)
158  {
159  motionState0_.t() = t;
160  motionState0_.deltaT() = deltaT;
161  }
162 
163  if (Pstream::master())
164  {
165  solver_->solve(tau, fx);
166  }
167 
168  Pstream::scatter(motionState_);
169 
170  // Update the body-state to correspond to the current joint-state
171  forwardDynamicsCorrection(motionState_);
172 }
173 
174 
175 void Foam::RBD::rigidBodyMotion::status(const label bodyID) const
176 {
177  const spatialTransform CofR(X0(bodyID));
178  const spatialVector vCofR(v(bodyID, Zero));
179 
180  Info<< "Rigid-body motion of the " << name(bodyID) << nl
181  << " Centre of rotation: " << CofR.r() << nl
182  << " Orientation: " << CofR.E() << nl
183  << " Linear velocity: " << vCofR.l() << nl
184  << " Angular velocity: " << vCofR.w()
185  << endl;
186 }
187 
188 
190 (
191  const label bodyID,
192  const pointField& initialPoints
193 ) const
194 {
195  // Calculate the transform from the initial state in the global frame
196  // to the current state in the global frame
197  spatialTransform X(X0(bodyID).inv() & X00(bodyID));
198 
199  tmp<pointField> tpoints(new pointField(initialPoints.size()));
200  pointField& points = tpoints.ref();
201 
202  forAll(points, i)
203  {
204  points[i] = X.transformPoint(initialPoints[i]);
205  }
206 
207  return tpoints;
208 }
209 
210 
212 (
213  const label bodyID,
214  const scalarField& weight,
215  const pointField& initialPoints
216 ) const
217 {
218  // Calculate the transform from the initial state in the global frame
219  // to the current state in the global frame
220  spatialTransform X(X0(bodyID).inv() & X00(bodyID));
221 
222  // Calculate the septernion equivalent of the transformation for 'slerp'
223  // interpolation
224  septernion s(X);
225 
226  tmp<pointField> tpoints(new pointField(initialPoints));
227  pointField& points = tpoints.ref();
228 
229  forAll(points, i)
230  {
231  // Move non-stationary points
232  if (weight[i] > SMALL)
233  {
234  // Use solid-body motion where weight = 1
235  if (weight[i] > 1 - SMALL)
236  {
237  points[i] = X.transformPoint(initialPoints[i]);
238  }
239  // Slerp septernion interpolation
240  else
241  {
242  points[i] =
243  slerp(septernion::I, s, weight[i])
244  .transformPoint(initialPoints[i]);
245  }
246  }
247  }
248 
249  return tpoints;
250 }
251 
252 
254 (
255  const labelList& bodyIDs,
256  const List<const scalarField*>& weights,
257  const pointField& initialPoints
258 ) const
259 {
260  List<septernion> ss(bodyIDs.size() + 1);
261  ss[bodyIDs.size()] = septernion::I;
262 
263  forAll(bodyIDs, bi)
264  {
265  const label bodyID = bodyIDs[bi];
266 
267  // Calculate the transform from the initial state in the global frame
268  // to the current state in the global frame
269  spatialTransform X(X0(bodyID).inv() & X00(bodyID));
270 
271  // Calculate the septernion equivalent of the transformation
272  ss[bi] = septernion(X);
273  }
274 
275  tmp<pointField> tpoints(new pointField(initialPoints));
276  pointField& points = tpoints.ref();
277 
278  List<scalar> w(ss.size());
279 
280  forAll(points, i)
281  {
282  // Initialize to 1 for the far-field weight
283  scalar sum1mw = 1;
284 
285  forAll(bodyIDs, bi)
286  {
287  w[bi] = (*(weights[bi]))[i];
288  sum1mw += w[bi]/(1 + SMALL - w[bi]);
289  }
290 
291  // Calculate the limiter for wi/(1 - wi) to ensure the sum(wi) = 1
292  scalar lambda = 1/sum1mw;
293 
294  // Limit wi/(1 - wi) and sum the resulting wi
295  scalar sumw = 0;
296  forAll(bodyIDs, bi)
297  {
298  w[bi] = lambda*w[bi]/(1 + SMALL - w[bi]);
299  sumw += w[bi];
300  }
301 
302  // Calculate the weight for the stationary far-field
303  w[bodyIDs.size()] = 1 - sumw;
304 
305  points[i] = average(ss, w).transformPoint(initialPoints[i]);
306  }
307 
308  return tpoints;
309 }
310 
311 
312 // ************************************************************************* //
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:147
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
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:175
Foam::UPstream::master
static bool master(const label communicator=worldComm)
Am I the master process.
Definition: UPstream.H:458
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:350
Foam::RBD::subBody::masterID
label masterID() const
Return the master body Id.
Definition: subBodyI.H:71
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:228
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 (uses stdout - output is on the master only)
Foam::name
word name(const complex &c)
Return string representation of complex.
Definition: complex.C:76
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
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:121
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:117
g
const uniformDimensionedVectorField & g
Definition: createFluidFields.H:24
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:385
Foam::Vector< scalar >
Foam::List< label >
points
const pointField & points
Definition: gmvOutputHeader.H:1
Foam::RBD::rigidBodyMotion::~rigidBodyMotion
~rigidBodyMotion()
Destructor.
Definition: rigidBodyMotion.C:110
Foam::septernion::I
static const septernion I
Definition: septernion.H:84
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:190
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:134
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