molecule.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2020 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 
69 (
70  moleculeCloud& cloud,
71  trackingData& td,
72  const scalar trackTime
73 )
74 {
75  td.switchProcessor = false;
76  td.keepParticle = true;
77 
78  const constantProperties& constProps(cloud.constProps(id_));
79 
81  {
82  // First leapfrog velocity adjust part, required before tracking+force
83  // part
84 
85  v_ += 0.5*trackTime*a_;
86 
87  pi_ += 0.5*trackTime*tau_;
88  }
89  else if (td.part() == trackingData::tpLinearTrack)
90  {
91  // Leapfrog tracking part
92 
93  while (td.keepParticle && !td.switchProcessor && stepFraction() < 1)
94  {
95  const scalar f = 1 - stepFraction();
96  trackToAndHitFace(f*trackTime*v_, f, cloud, td);
97  }
98  }
99  else if (td.part() == trackingData::tpRotationalTrack)
100  {
101  // Leapfrog orientation adjustment, carried out before force calculation
102  // but after tracking stage, i.e. rotation carried once linear motion
103  // complete.
104 
105  if (!constProps.pointMolecule())
106  {
107  const diagTensor& momentOfInertia(constProps.momentOfInertia());
108 
109  tensor R;
110 
111  if (!constProps.linearMolecule())
112  {
113  R = rotationTensorX(0.5*trackTime*pi_.x()/momentOfInertia.xx());
114  pi_ = pi_ & R;
115  Q_ = Q_ & R;
116  }
117 
118  R = rotationTensorY(0.5*trackTime*pi_.y()/momentOfInertia.yy());
119  pi_ = pi_ & R;
120  Q_ = Q_ & R;
121 
122  R = rotationTensorZ(trackTime*pi_.z()/momentOfInertia.zz());
123  pi_ = pi_ & R;
124  Q_ = Q_ & R;
125 
126  R = rotationTensorY(0.5*trackTime*pi_.y()/momentOfInertia.yy());
127  pi_ = pi_ & R;
128  Q_ = Q_ & R;
129 
130  if (!constProps.linearMolecule())
131  {
132  R = rotationTensorX(0.5*trackTime*pi_.x()/momentOfInertia.xx());
133  pi_ = pi_ & R;
134  Q_ = Q_ & R;
135  }
136  }
137 
138  setSitePositions(constProps);
139  }
140  else if (td.part() == trackingData::tpVelocityHalfStep1)
141  {
142  // Second leapfrog velocity adjust part, required after tracking+force
143  // part
144 
145  scalar m = constProps.mass();
146 
147  a_ = Zero;
148 
149  tau_ = Zero;
150 
151  forAll(siteForces_, s)
152  {
153  const vector& f = siteForces_[s];
154 
155  a_ += f/m;
156 
157  tau_ += (constProps.siteReferencePositions()[s] ^ (Q_.T() & f));
158  }
159 
160  v_ += 0.5*trackTime*a_;
161 
162  pi_ += 0.5*trackTime*tau_;
163 
164  if (constProps.pointMolecule())
165  {
166  tau_ = Zero;
167 
168  pi_ = Zero;
169  }
170 
171  if (constProps.linearMolecule())
172  {
173  tau_.x() = 0.0;
174 
175  pi_.x() = 0.0;
176  }
177  }
178  else
179  {
181  << td.part() << " is an invalid part of the integration method."
182  << abort(FatalError);
183  }
184 
185  return td.keepParticle;
186 }
187 
188 
190 {
192 
193  Q_ = transform.T() & Q_;
194 
195  v_ = transform.transform(v_);
196 
197  a_ = transform.transform(a_);
198 
199  pi_ = Q_.T() & transform.transform(Q_ & pi_);
200 
201  tau_ = Q_.T() & transform.transform(Q_ & tau_);
202 
203  rf_ = transform.transform(rf_);
204 
205  transform.transformList(siteForces_);
206 
207  sitePositions_ = transform.transformPosition(vectorField(sitePositions_));
208 
209  if (special_ == SPECIAL_TETHERED)
210  {
211  specialPosition_ = transform.transformPosition(specialPosition_);
212  }
213 }
214 
215 
217 {
218  sitePositions_ = position() + (Q_ & constProps.siteReferencePositions());
219 }
220 
221 
223 {
224  sitePositions_.setSize(size);
225 
226  siteForces_.setSize(size);
227 }
228 
229 
231 {
232  return false;
233 }
234 
235 
237 {
238  td.switchProcessor = true;
239 }
240 
241 
243 {
244  const vector nw = normal();
245  const scalar vn = v_ & nw;
246 
247  // Specular reflection
248  if (vn > 0)
249  {
250  v_ -= 2*vn*nw;
251  }
252 }
253 
254 
255 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
Vector-tensor class used to perform translations, rotations and scaling operations in 3D space...
Definition: transformer.H:83
Class to hold molecule constant properties.
Definition: molecule.H:90
void setSitePositions(const constantProperties &constProps)
Definition: molecule.C:216
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
bool hitPatch(moleculeCloud &cloud, trackingData &td)
Overridable function to handle the particle hitting a patch.
Definition: molecule.C:230
Tensor< Cmpt > T() const
Return transpose.
Definition: TensorI.H:331
scalar trackToAndHitFace(const vector &displacement, const scalar fraction, TrackCloudType &cloud, trackingData &td)
Convenience function. Combines trackToFace and hitFace.
void hitProcessorPatch(moleculeCloud &cloud, trackingData &td)
Overridable function to handle the particle hitting a processorPatch.
Definition: molecule.C:236
virtual void transformProperties(const transformer &)
Transform the physical properties of the particle.
Definition: molecule.C:189
virtual void transformProperties(const transformer &)
Transform the physical properties of the particle.
Definition: particle.C:1048
const Cmpt & y() const
Definition: VectorI.H:81
bool move(moleculeCloud &, trackingData &, const scalar trackTime)
Definition: molecule.C:69
vector normal() const
Return the normal of the tri on tetFacei_ for the.
Definition: particleI.H:248
void transformList(Container< Type > &) const
Transform the given container.
const Field< vector > & siteReferencePositions() const
Definition: moleculeI.H:366
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))
dimensionedScalar cos(const dimensionedScalar &ds)
bool switchProcessor
Flag to switch processor.
Definition: particle.H:109
trackPart part() const
Return the part of the tracking operation taking place.
Definition: molecule.H:203
static const zero Zero
Definition: zero.H:97
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const Cmpt & x() const
Definition: VectorI.H:75
void hitWallPatch(moleculeCloud &cloud, trackingData &td)
Overridable function to handle the particle hitting a wallPatch.
Definition: molecule.C:242
dimensionedScalar sin(const dimensionedScalar &ds)
const tensor & T() const
Return the transformation tensor.
Definition: transformerI.H:94
labelList f(nPoints)
scalar stepFraction() const
Return the fraction of time-step completed.
Definition: particleI.H:161
void setSize(const label)
Reset size of List.
Definition: List.C:281
#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
Class used to pass tracking data to the trackToFace function.
Definition: molecule.H:166
void setSiteSizes(label size)
Definition: molecule.C:222
Field< vector > vectorField
Specialisation of Field<T> for vector.
Type transform(const Type &) const
Transform the given type.
const List< molecule::constantProperties > constProps() const
vector transformPosition(const vector &v) const
Transform the given position.
Definition: transformerI.H:153
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.
vector position() const
Return current particle position.
Definition: particleI.H:278
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
label nw
Definition: createFields.H:12
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:477