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-------------------------------------------------------------------------------
11License
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
35void 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
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
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 (
71 Function1<scalar>::NewIfPresent
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
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 (
105 Function1<scalar>::NewIfPresent
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
200void 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
214const 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
222const 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// ************************************************************************* //
Top level data entry class for use in dictionaries. Provides a mechanism to specify a variable as a c...
Definition: Function1.H:96
static void scatter(const List< commsStruct > &comms, T &value, const int tag, const label comm)
Broadcast data: Distribute without modification.
Holds the motion state of rigid-body model.
const scalarField & qDdot() const
Return access to the joint acceleration.
Basic rigid-body model representing a system of rigid-bodies connected by 1-6 DoF joints.
DynamicList< spatialTransform > X0_
Transform for external forces to the bodies reference frame.
void forwardDynamicsCorrection(const rigidBodyModelState &state) const
Correct the velocity and acceleration of the bodies in the model.
const vector & g() const
Return the acceleration due to gravity.
void forwardDynamics(rigidBodyModelState &state, const scalarField &tau, const Field< spatialVector > &fx) const
Calculate the joint acceleration qDdot from the joint state q,.
Six degree of freedom motion for a rigid body.
spatialTransform X00(const label bodyId) const
Return the initial transform to the global frame for the.
tmp< pointField > transformPoints(const label bodyID, const pointField &initialPoints) const
Transform the given initial pointField of the specified body.
const vector vCofR(const label bodyID) const
Report linear velocity of the given body.
void forwardDynamics(rigidBodyModelState &state, const scalarField &tau, const Field< spatialVector > &fx) const
Calculate and optionally relax the joint acceleration qDdot from.
const vector cCofR(const label bodyID) const
Report CofR of the given body.
const spatialTransform & masterXT() const
Return the transform with respect to the master body.
Definition: subBodyI.H:77
label masterID() const
Return the master body Id.
Definition: subBodyI.H:71
Vector< Cmpt > w() const
Return the angular part of the spatial vector as a vector.
Vector< Cmpt > l() const
Return the linear part of the spatial vector as a vector.
A simple wrapper around bool so that it can be read as a word: true/false, on/off,...
Definition: Switch.H:78
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:80
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
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
bool status() const
Return the execution status (on/off) of the function objects.
Septernion class used to perform translations and rotations in 3D space.
Definition: septernion.H:67
static const septernion I
Definition: septernion.H:84
Compact representation of the Plücker spatial transformation tensor in terms of the rotation tensor E...
const vector & r() const
Return the translation vector.
const tensor & E() const
Return the rotation tensor.
vector transformPoint(const vector &p) const
Transform position p.
splitCell * master() const
Definition: splitCell.H:113
A class for managing temporary objects.
Definition: tmp.H:65
T & ref() const
Definition: tmpI.H:227
A class for handling words, derived from Foam::string.
Definition: word.H:68
const pointField & points
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))
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh > > &tdf1, const word &name, const dimensionSet &dimensions)
Global function forwards to reuseTmpDimensionedField::New.
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
quaternion slerp(const quaternion &qa, const quaternion &qb, const scalar t)
Spherical linear interpolation of quaternions.
Definition: quaternion.C:82
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
scalarList X0(nSpecie, Zero)
dictionary dict
CEqn solve()
dimensionedScalar lambda("lambda", dimTime/sqr(dimLength), laminarTransport)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333