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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2017 OpenFOAM Foundation
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "moleculeCloud.H"
29 #include "molecule.H"
30 #include "Random.H"
31 #include "Time.H"
32 
33 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
34 
35 Foam::tensor Foam::molecule::rotationTensorX(scalar phi) const
36 {
37  return tensor
38  (
39  1, 0, 0,
40  0, Foam::cos(phi), -Foam::sin(phi),
42  );
43 }
44 
45 
46 Foam::tensor Foam::molecule::rotationTensorY(scalar phi) const
47 {
48  return tensor
49  (
50  Foam::cos(phi), 0, Foam::sin(phi),
51  0, 1, 0,
52  -Foam::sin(phi), 0, Foam::cos(phi)
53  );
54 }
55 
56 
57 Foam::tensor Foam::molecule::rotationTensorZ(scalar phi) const
58 {
59  return tensor
60  (
61  Foam::cos(phi), -Foam::sin(phi), 0,
62  Foam::sin(phi), Foam::cos(phi), 0,
63  0, 0, 1
64  );
65 }
66 
67 
68 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
69 
71 (
73  trackingData& td,
74  const scalar trackTime
75 )
76 {
77  td.switchProcessor = false;
78  td.keepParticle = true;
79 
80  const constantProperties& constProps(cloud.constProps(id_));
81 
82  if (td.part() == 0)
83  {
84  // First leapfrog velocity adjust part, required before tracking+force
85  // part
86 
87  v_ += 0.5*trackTime*a_;
88 
89  pi_ += 0.5*trackTime*tau_;
90  }
91  else if (td.part() == 1)
92  {
93  // Leapfrog tracking part
94 
95  while (td.keepParticle && !td.switchProcessor && stepFraction() < 1)
96  {
97  const scalar f = 1 - stepFraction();
98  trackToAndHitFace(f*trackTime*v_, f, cloud, td);
99  }
100  }
101  else if (td.part() == 2)
102  {
103  // Leapfrog orientation adjustment, carried out before force calculation
104  // but after tracking stage, i.e. rotation carried once linear motion
105  // complete.
106 
107  if (!constProps.pointMolecule())
108  {
109  const diagTensor& momentOfInertia(constProps.momentOfInertia());
110 
111  tensor R;
112 
113  if (!constProps.linearMolecule())
114  {
115  R = rotationTensorX(0.5*trackTime*pi_.x()/momentOfInertia.xx());
116  pi_ = pi_ & R;
117  Q_ = Q_ & R;
118  }
119 
120  R = rotationTensorY(0.5*trackTime*pi_.y()/momentOfInertia.yy());
121  pi_ = pi_ & R;
122  Q_ = Q_ & R;
123 
124  R = rotationTensorZ(trackTime*pi_.z()/momentOfInertia.zz());
125  pi_ = pi_ & R;
126  Q_ = Q_ & R;
127 
128  R = rotationTensorY(0.5*trackTime*pi_.y()/momentOfInertia.yy());
129  pi_ = pi_ & R;
130  Q_ = Q_ & R;
131 
132  if (!constProps.linearMolecule())
133  {
134  R = rotationTensorX(0.5*trackTime*pi_.x()/momentOfInertia.xx());
135  pi_ = pi_ & R;
136  Q_ = Q_ & R;
137  }
138  }
139 
140  setSitePositions(constProps);
141  }
142  else if (td.part() == 3)
143  {
144  // Second leapfrog velocity adjust part, required after tracking+force
145  // part
146 
147  scalar m = constProps.mass();
148 
149  a_ = Zero;
150 
151  tau_ = Zero;
152 
153  forAll(siteForces_, s)
154  {
155  const vector& f = siteForces_[s];
156 
157  a_ += f/m;
158 
159  tau_ += (constProps.siteReferencePositions()[s] ^ (Q_.T() & f));
160  }
161 
162  v_ += 0.5*trackTime*a_;
163 
164  pi_ += 0.5*trackTime*tau_;
165 
166  if (constProps.pointMolecule())
167  {
168  tau_ = Zero;
169 
170  pi_ = Zero;
171  }
172 
173  if (constProps.linearMolecule())
174  {
175  tau_.x() = 0.0;
176 
177  pi_.x() = 0.0;
178  }
179  }
180  else
181  {
183  << td.part() << " is an invalid part of the integration method."
184  << abort(FatalError);
185  }
186 
187  return td.keepParticle;
188 }
189 
190 
192 {
194 
195  Q_ = T & Q_;
196 
197  v_ = transform(T, v_);
198 
199  a_ = transform(T, a_);
200 
201  pi_ = Q_.T() & transform(T, Q_ & pi_);
202 
203  tau_ = Q_.T() & transform(T, Q_ & tau_);
204 
205  rf_ = transform(T, rf_);
206 
207  sitePositions_ = position() + (T & (sitePositions_ - position()));
208 
209  siteForces_ = T & siteForces_;
210 }
211 
212 
214 {
215  particle::transformProperties(separation);
216 
217  if (special_ == SPECIAL_TETHERED)
218  {
219  specialPosition_ += separation;
220  }
221 
222  sitePositions_ = sitePositions_ + separation;
223 }
224 
225 
227 {
228  sitePositions_ = position() + (Q_ & constProps.siteReferencePositions());
229 }
230 
231 
233 {
234  sitePositions_.setSize(size);
235 
236  siteForces_.setSize(size);
237 }
238 
239 
241 {
242  return false;
243 }
244 
245 
247 {
248  td.switchProcessor = true;
249 }
250 
251 
253 {
254  const vector nw = normal();
255 
256  const scalar vn = v_ & nw;
257 
258  // Specular reflection
259  if (vn > 0)
260  {
261  v_ -= 2*vn*nw;
262  }
263 }
264 
265 
266 // ************************************************************************* //
Foam::Tensor< scalar >
Foam::DiagTensor< scalar >
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::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::sin
dimensionedScalar sin(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:264
Foam::transform
dimensionSet transform(const dimensionSet &ds)
Return the argument; transformations do not change the dimensions.
Definition: dimensionSet.C:521
Foam::molecule::trackingData::part
label part() const
Definition: molecule.H:187
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::molecule::hitWallPatch
void hitWallPatch(moleculeCloud &cloud, trackingData &td)
Overridable function to handle the particle hitting a wallPatch.
Definition: molecule.C:252
Foam::molecule::setSitePositions
void setSitePositions(const constantProperties &constProps)
Definition: molecule.C:226
moleculeCloud.H
R
#define R(A, B, C, D, E, F, K, M)
Foam::molecule::constantProperties::siteReferencePositions
const Field< vector > & siteReferencePositions() const
Definition: moleculeI.H:374
Foam::T
void T(FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
Definition: FieldFieldFunctions.C:58
Foam::molecule::trackingData
Class used to pass tracking data to the trackToFace function.
Definition: molecule.H:167
Foam::particle::trackingData::keepParticle
bool keepParticle
Flag to indicate whether to keep particle (false = delete)
Definition: particle.H:105
Foam::particle::trackingData::switchProcessor
bool switchProcessor
Flag to switch processor.
Definition: particle.H:102
phi
surfaceScalarField & phi
Definition: setRegionFluidFields.H:8
Foam::molecule::setSiteSizes
void setSiteSizes(label size)
Definition: molecule.C:232
Foam::molecule::hitProcessorPatch
void hitProcessorPatch(moleculeCloud &cloud, trackingData &td)
Overridable function to handle the particle hitting a processorPatch.
Definition: molecule.C:246
Foam::FatalError
error FatalError
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:144
Random.H
Foam::molecule::move
bool move(moleculeCloud &, trackingData &, const scalar trackTime)
Definition: molecule.C:71
Time.H
Foam::cloud
A cloud is a registry collection of lagrangian particles.
Definition: cloud.H:57
Foam::momentOfInertia
Calculates the inertia tensor and principal axes and moments of a polyhedra/cells/triSurfaces....
Definition: momentOfInertia.H:65
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
f
labelList f(nPoints)
Foam::moleculeCloud
Definition: moleculeCloud.H:58
Foam::particle::transformProperties
virtual void transformProperties(const tensor &T)
Transform the physical properties of the particle.
Definition: particle.C:1072
Foam::Vector< scalar >
Foam::molecule::constantProperties
Class to hold molecule constant properties.
Definition: molecule.H:91
molecule.H
Foam::molecule::hitPatch
bool hitPatch(moleculeCloud &cloud, trackingData &td)
Overridable function to handle the particle hitting a patch.
Definition: molecule.C:240
Foam::tensor
Tensor< scalar > tensor
Tensor of scalars, i.e. Tensor<scalar>.
Definition: symmTensor.H:61
Foam::molecule::transformProperties
virtual void transformProperties(const tensor &T)
Transform the physical properties of the particle.
Definition: molecule.C:191
Foam::cos
dimensionedScalar cos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:265