molecule.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 | Copyright (C) 2011-2013 OpenFOAM Foundation
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 \*---------------------------------------------------------------------------*/
25 
26 #include "moleculeCloud.H"
27 #include "molecule.H"
28 #include "Random.H"
29 #include "Time.H"
30 
31 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
32 
33 Foam::tensor Foam::molecule::rotationTensorX(scalar phi) const
34 {
35  return tensor
36  (
37  1, 0, 0,
38  0, Foam::cos(phi), -Foam::sin(phi),
39  0, Foam::sin(phi), Foam::cos(phi)
40  );
41 }
42 
43 
44 Foam::tensor Foam::molecule::rotationTensorY(scalar phi) const
45 {
46  return tensor
47  (
48  Foam::cos(phi), 0, Foam::sin(phi),
49  0, 1, 0,
50  -Foam::sin(phi), 0, Foam::cos(phi)
51  );
52 }
53 
54 
55 Foam::tensor Foam::molecule::rotationTensorZ(scalar phi) const
56 {
57  return tensor
58  (
59  Foam::cos(phi), -Foam::sin(phi), 0,
60  Foam::sin(phi), Foam::cos(phi), 0,
61  0, 0, 1
62  );
63 }
64 
65 
66 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
67 
68 bool Foam::molecule::move(molecule::trackingData& td, const scalar trackTime)
69 {
70  td.switchProcessor = false;
71  td.keepParticle = true;
72 
73  const constantProperties& constProps(td.cloud().constProps(id_));
74 
75  if (td.part() == 0)
76  {
77  // First leapfrog velocity adjust part, required before tracking+force
78  // part
79 
80  v_ += 0.5*trackTime*a_;
81 
82  pi_ += 0.5*trackTime*tau_;
83  }
84  else if (td.part() == 1)
85  {
86  // Leapfrog tracking part
87 
88  scalar tEnd = (1.0 - stepFraction())*trackTime;
89  scalar dtMax = tEnd;
90 
91  while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
92  {
93  // set the lagrangian time-step
94  scalar dt = min(dtMax, tEnd);
95 
96  dt *= trackToFace(position() + dt*v_, td);
97 
98  tEnd -= dt;
99  stepFraction() = 1.0 - tEnd/trackTime;
100  }
101  }
102  else if (td.part() == 2)
103  {
104  // Leapfrog orientation adjustment, carried out before force calculation
105  // but after tracking stage, i.e. rotation carried once linear motion
106  // complete.
107 
108  if (!constProps.pointMolecule())
109  {
110  const diagTensor& momentOfInertia(constProps.momentOfInertia());
111 
112  tensor R;
113 
114  if (!constProps.linearMolecule())
115  {
116  R = rotationTensorX(0.5*trackTime*pi_.x()/momentOfInertia.xx());
117  pi_ = pi_ & R;
118  Q_ = Q_ & R;
119  }
120 
121  R = rotationTensorY(0.5*trackTime*pi_.y()/momentOfInertia.yy());
122  pi_ = pi_ & R;
123  Q_ = Q_ & R;
124 
125  R = rotationTensorZ(trackTime*pi_.z()/momentOfInertia.zz());
126  pi_ = pi_ & R;
127  Q_ = Q_ & R;
128 
129  R = rotationTensorY(0.5*trackTime*pi_.y()/momentOfInertia.yy());
130  pi_ = pi_ & R;
131  Q_ = Q_ & R;
132 
133  if (!constProps.linearMolecule())
134  {
135  R = rotationTensorX(0.5*trackTime*pi_.x()/momentOfInertia.xx());
136  pi_ = pi_ & R;
137  Q_ = Q_ & R;
138  }
139  }
140 
141  setSitePositions(constProps);
142  }
143  else if (td.part() == 3)
144  {
145  // Second leapfrog velocity adjust part, required after tracking+force
146  // part
147 
148  scalar m = constProps.mass();
149 
150  a_ = vector::zero;
151 
152  tau_ = vector::zero;
153 
154  forAll(siteForces_, s)
155  {
156  const vector& f = siteForces_[s];
157 
158  a_ += f/m;
159 
160  tau_ += (constProps.siteReferencePositions()[s] ^ (Q_.T() & f));
161  }
162 
163  v_ += 0.5*trackTime*a_;
164 
165  pi_ += 0.5*trackTime*tau_;
166 
167  if (constProps.pointMolecule())
168  {
169  tau_ = vector::zero;
170 
171  pi_ = vector::zero;
172  }
173 
174  if (constProps.linearMolecule())
175  {
176  tau_.x() = 0.0;
177 
178  pi_.x() = 0.0;
179  }
180  }
181  else
182  {
183  FatalErrorIn("molecule::move(trackingData&, const scalar)") << nl
184  << td.part() << " is an invalid part of the integration method."
185  << abort(FatalError);
186  }
187 
188  return td.keepParticle;
189 }
190 
191 
193 {
195 
196  Q_ = T & Q_;
197 
198  v_ = transform(T, v_);
199 
200  a_ = transform(T, a_);
201 
202  pi_ = Q_.T() & transform(T, Q_ & pi_);
203 
204  tau_ = Q_.T() & transform(T, Q_ & tau_);
205 
206  rf_ = transform(T, rf_);
207 
208  sitePositions_ = position_ + (T & (sitePositions_ - position_));
209 
210  siteForces_ = T & siteForces_;
211 }
212 
213 
215 {
216  particle::transformProperties(separation);
217 
218  if (special_ == SPECIAL_TETHERED)
219  {
220  specialPosition_ += separation;
221  }
222 
223  sitePositions_ = sitePositions_ + separation;
224 }
225 
226 
228 {
229  sitePositions_ = position_ + (Q_ & constProps.siteReferencePositions());
230 }
231 
232 
234 {
235  sitePositions_.setSize(size);
236 
237  siteForces_.setSize(size);
238 }
239 
240 
242 (
243  const polyPatch&,
244  trackingData&,
245  const label,
246  const scalar,
247  const tetIndices&
248 )
249 {
250  return false;
251 }
252 
253 
255 (
256  const processorPolyPatch&,
257  trackingData& td
258 )
259 {
260  td.switchProcessor = true;
261 }
262 
263 
265 (
266  const wallPolyPatch& wpp,
267  trackingData& td,
268  const tetIndices& tetIs
269 )
270 {
271  // Use of the normal from tetIs is not required as
272  // hasWallImpactDistance for a moleculeCloud is false.
273  vector nw = normal();
274  nw /= mag(nw);
275 
276  scalar vn = v_ & nw;
277 
278  // Specular reflection
279  if (vn > 0)
280  {
281  v_ -= 2*vn*nw;
282  }
283 }
284 
285 
287 (
288  const polyPatch&,
289  trackingData& td
290 )
291 {
292  td.keepParticle = false;
293 }
294 
295 
296 // ************************************************************************* //
const Field< vector > & siteReferencePositions() const
Definition: moleculeI.H:331
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().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 ))
dimensioned< scalar > mag(const dimensioned< Type > &)
labelList f(nPoints)
void setSitePositions(const constantProperties &constProps)
Definition: molecule.C:227
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
#define R(A, B, C, D, E, F, K, M)
bool keepParticle
Flag to indicate whether to keep particle (false = delete)
Definition: particle.H:112
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Class used to pass tracking data to the trackToFace function.
Definition: molecule.H:157
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
Class to hold molecule constant properties.
Definition: molecule.H:81
label nw
Definition: createFields.H:25
void hitProcessorPatch(const processorPolyPatch &, trackingData &td)
Overridable function to handle the particle hitting a processorPatch.
Definition: molecule.C:255
bool move(trackingData &, const scalar trackTime)
Definition: molecule.C:68
static const char nl
Definition: Ostream.H:260
void setSize(const label)
Reset size of List.
Definition: List.C:318
const Cmpt & y() const
Definition: VectorI.H:71
Neighbour processor patch.
vector normal() const
Return the normal of the tri on tetFaceI_ for the.
Definition: particleI.H:663
const vector & position() const
Return current particle position.
Definition: particleI.H:603
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:465
#define forAll(list, i)
Definition: UList.H:421
dimensionedScalar cos(const dimensionedScalar &ds)
const Cmpt & x() const
Definition: VectorI.H:65
Calculates the inertia tensor and principal axes and moments of a polyhedra/cells/triSurfaces. Inertia can either be of the solid body or of a thin shell.
bool hitPatch(const polyPatch &, trackingData &td, const label patchI, const scalar trackFraction, const tetIndices &tetIs)
Overridable function to handle the particle hitting a patch.
Definition: molecule.C:242
scalar trackToFace(const vector &endPosition, TrackData &td)
Track particle to a given position and returns 1.0 if the.
CloudType & cloud()
Return a reference to the cloud.
Definition: particle.H:125
errorManip< error > abort(error &err)
Definition: errorManip.H:131
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
scalar & stepFraction()
Return the fraction of time-step completed.
Definition: particleI.H:834
virtual void transformProperties(const tensor &T)
Transform the physical properties of the particle.
Definition: molecule.C:192
void hitWallPatch(const wallPolyPatch &, trackingData &td, const tetIndices &)
Overridable function to handle the particle hitting a wallPatch.
Definition: molecule.C:265
Foam::wallPolyPatch.
Definition: wallPolyPatch.H:48
bool switchProcessor
Flag to switch processor.
Definition: particle.H:109
error FatalError
const List< molecule::constantProperties > constProps() const
static const Vector zero
Definition: Vector.H:80
Tensor< Cmpt > T() const
Transpose.
Definition: TensorI.H:286
void setSiteSizes(label size)
Definition: molecule.C:233
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:73
virtual void transformProperties(const tensor &T)
Transform the physical properties of the particle.
Definition: particle.C:123
vector position_
Position of particle.
Definition: particle.H:140
dimensionedScalar sin(const dimensionedScalar &ds)