FreeFOAM The Cross-Platform CFD Toolkit
sixDoFRigidBodyMotion.H
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2009-2010 OpenCFD Ltd.
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 Class
25  Foam::sixDoFRigidBodyMotion
26 
27 Description
28  Six degree of freedom motion for a rigid body. Angular momentum
29  stored in body fixed reference frame. Reference orientation of
30  the body (where Q = I) must align with the cartesian axes such
31  that the Inertia tensor is in principle component form.
32 
33  Symplectic motion as per:
34 
35  title = {Symplectic splitting methods for rigid body molecular dynamics},
36  publisher = {AIP},
37  year = {1997},
38  journal = {The Journal of Chemical Physics},
39  volume = {107},
40  number = {15},
41  pages = {5840-5851},
42  url = {http://link.aip.org/link/?JCP/107/5840/1},
43  doi = {10.1063/1.474310}
44 
45  Can add restraints (i.e. a spring) and constraints (i.e. motion
46  may only be on a plane).
47 
48 SourceFiles
49  sixDoFRigidBodyMotionI.H
50  sixDoFRigidBodyMotion.C
51  sixDoFRigidBodyMotionIO.C
52 
53 \*---------------------------------------------------------------------------*/
54 
55 #ifndef sixDoFRigidBodyMotion_H
56 #define sixDoFRigidBodyMotion_H
57 
59 #include <OpenFOAM/pointField.H>
62 
63 
64 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
65 
66 namespace Foam
67 {
68 
69 // Forward declaration of classes
70 class Istream;
71 class Ostream;
72 
73 // Forward declaration of friend functions and operators
74 class sixDoFRigidBodyMotion;
75 Istream& operator>>(Istream&, sixDoFRigidBodyMotion&);
76 Ostream& operator<<(Ostream&, const sixDoFRigidBodyMotion&);
77 
78 
79 /*---------------------------------------------------------------------------*\
80  Class sixDoFRigidBodyMotion Declaration
81 \*---------------------------------------------------------------------------*/
82 
84 {
85  // Private data
86 
87  //- Motion state data object
88  sixDoFRigidBodyMotionState motionState_;
89 
90  //- Restraints on the motion
92 
93  //- Names of the restraints
94  wordList restraintNames_;
95 
96  //- Constaints on the motion
98 
99  //- Names of the constraints
100  wordList constraintNames_;
101 
102  //- Maximum number of iterations allowed to attempt to obey
103  // constraints
104  label maxConstraintIterations_;
105 
106  //- Centre of mass of initial state
107  point initialCentreOfMass_;
108 
109  //- Orientation of initial state
110  tensor initialQ_;
111 
112  //- Moment of inertia of the body in reference configuration
113  // (Q = I)
114  diagTensor momentOfInertia_;
115 
116  //- Mass of the body
117  scalar mass_;
118 
119  //- Acceleration damping coefficient. Modify applied acceleration:
120  // v1 = v0 + a*dt - cDamp*a*dt
121  // = v0 + dt*f*(1 - cDamp)/m
122  // Increases effective mass by 1/(1 - cDamp).
123  scalar cDamp_;
124 
125  //- Acceleration magnitude limit - clips large accelerations
126  scalar aLim_;
127 
128  //- Switch to turn reporting of motion data on and off
129  Switch report_;
130 
131 
132  // Private Member Functions
133 
134  //- Calculate the rotation tensor around the body reference
135  // frame x-axis by the given angle
136  inline tensor rotationTensorX(scalar deltaT) const;
137 
138  //- Calculate the rotation tensor around the body reference
139  // frame y-axis by the given angle
140  inline tensor rotationTensorY(scalar deltaT) const;
141 
142  //- Calculate the rotation tensor around the body reference
143  // frame z-axis by the given angle
144  inline tensor rotationTensorZ(scalar deltaT) const;
145 
146  //- Apply rotation tensors to Q for the given torque (pi) and deltaT
147  inline void rotate(tensor& Q, vector& pi, scalar deltaT) const;
148 
149  //- Apply the restraints to the object
150  void applyRestraints();
151 
152  //- Apply the constraints to the object
153  void applyConstraints(scalar deltaT);
154 
155  // Access functions retained as private because of the risk of
156  // confusion over what is a body local frame vector and what is global
157 
158  // Access
159 
160  //- Return access to the motion state
161  inline const sixDoFRigidBodyMotionState& motionState() const;
162 
163  //- Return access to the restraints
165  restraints() const;
166 
167  //- Return access to the restraintNames
168  inline const wordList& restraintNames() const;
169 
170  //- Return access to the constraints
172  constraints() const;
173 
174  //- Return access to the constraintNames
175  inline const wordList& constraintNames() const;
176 
177  //- Return access to the maximum allowed number of
178  // constraint iterations
179  inline label maxConstraintIterations() const;
180 
181  //- Return access to the initial centre of mass
182  inline const point& initialCentreOfMass() const;
183 
184  //- Return access to the initial orientation
185  inline const tensor& initialQ() const;
186 
187  //- Return access to the orientation
188  inline const tensor& Q() const;
189 
190  //- Return access to velocity
191  inline const vector& v() const;
192 
193  //- Return access to acceleration
194  inline const vector& a() const;
195 
196  //- Return access to angular momentum
197  inline const vector& pi() const;
198 
199  //- Return access to torque
200  inline const vector& tau() const;
201 
202 
203  // Edit
204 
205  //- Return access to the centre of mass
206  inline point& initialCentreOfMass();
207 
208  //- Return access to the centre of mass
209  inline tensor& initialQ();
210 
211  //- Return non-const access to the orientation
212  inline tensor& Q();
213 
214  //- Return non-const access to vector
215  inline vector& v();
216 
217  //- Return non-const access to acceleration
218  inline vector& a();
219 
220  //- Return non-const access to angular momentum
221  inline vector& pi();
222 
223  //- Return non-const access to torque
224  inline vector& tau();
225 
226 
227 public:
228 
229  // Constructors
230 
231  //- Construct null
233 
234  //- Construct from components
236  (
237  const point& centreOfMass,
238  const tensor& Q,
239  const vector& v,
240  const vector& a,
241  const vector& pi,
242  const vector& tau,
243  scalar mass,
244  const point& initialCentreOfMass,
245  const tensor& initialQ,
247  scalar cDamp = 0.0,
248  scalar aLim = VGREAT,
249  bool report = false
250  );
251 
252  //- Construct from dictionary
253  sixDoFRigidBodyMotion(const dictionary& dict);
254 
255  //- Construct as copy
257 
258 
259  //- Destructor
261 
262 
263  // Member Functions
264 
265  //- Add restraints to the motion, public to allow external
266  // addition of restraints after construction
267  void addRestraints(const dictionary& dict);
268 
269  //- Add restraints to the motion, public to allow external
270  // addition of restraints after construction
271  void addConstraints(const dictionary& dict);
272 
273  //- First leapfrog velocity adjust and motion part, required
274  // before force calculation. Takes old timestep for variable
275  // timestep cases.
276  void updatePosition
277  (
278  scalar deltaT,
279  scalar deltaT0
280  );
281 
282  //- Second leapfrog velocity adjust part, required after motion and
283  // force calculation
284  void updateForce
285  (
286  const vector& fGlobal,
287  const vector& tauGlobal,
288  scalar deltaT
289  );
290 
291  //- Global forces supplied at locations, calculating net force
292  // and moment
293  void updateForce
294  (
295  const pointField& positions,
296  const vectorField& forces,
297  scalar deltaT
298  );
299 
300  //- Transform the given initial state pointField by the current
301  // motion state
303  (
304  const pointField& pInitial
305  ) const;
306 
307  //- Transform the given initial state point by the current motion
308  // state
309  inline point currentPosition(const point& pInitial) const;
310 
311  //- Transform the given initial state direction by the current
312  // motion state
313  inline vector currentOrientation(const vector& vInitial) const;
314 
315  //- Access the orientation tensor, Q.
316  // globalVector = Q & bodyLocalVector
317  // bodyLocalVector = Q.T() & globalVector
318  inline const tensor& orientation() const;
319 
320  //- Predict the position of the supplied initial state point
321  // after deltaT given the current motion state and the
322  // additional supplied force and moment
324  (
325  const point& pInitial,
326  const vector& deltaForce,
327  const vector& deltaMoment,
328  scalar deltaT
329  ) const;
330 
331  //- Predict the orientation of the supplied initial state
332  // vector after deltaT given the current motion state and the
333  // additional supplied moment
335  (
336  const vector& vInitial,
337  const vector& deltaMoment,
338  scalar deltaT
339  ) const;
340 
341  //- Return the angular velocity in the global frame
342  inline vector omega() const;
343 
344  //- Return the velocity of a position given by the current
345  // motion state
346  inline point currentVelocity(const point& pt) const;
347 
348  //- Report the status of the motion
349  void status() const;
350 
351 
352  // Access
353 
354  //- Return const access to the centre of mass
355  inline const point& centreOfMass() const;
356 
357  //- Return access to the inertia tensor
358  inline const diagTensor& momentOfInertia() const;
359 
360  //- Return const access to the mass
361  inline scalar mass() const;
362 
363  //- Return the report Switch
364  inline bool report() const;
365 
366 
367  // Edit
368 
369  //- Return non-const access to the centre of mass
370  inline point& centreOfMass();
371 
372  //- Return non-const access to the inertia tensor
373  inline diagTensor& momentOfInertia();
374 
375  //- Return non-const access to the mass
376  inline scalar& mass();
377 
378 
379  //- Write
380  void write(Ostream&) const;
381 
382 
383  // IOstream Operators
384 
387 };
388 
389 
390 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
391 
392 } // End namespace Foam
393 
394 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
395 
396 #include "sixDoFRigidBodyMotionI.H"
397 
398 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
399 
400 #endif
401 
402 // ************************ vim: set sw=4 sts=4 et: ************************ //