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-2024 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 "randomGenerator.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 (
71  trackingData& td
72 )
73 {
74  td.keepParticle = true;
75  td.sendToProc = -1;
76 
77  const constantProperties& constProps(cloud.constProps(id_));
78 
79  const scalar trackTime = td.mesh.time().deltaTValue();
80 
81  if (td.part() == trackingData::tpVelocityHalfStep0)
82  {
83  // First leapfrog velocity adjust part, required before tracking+force
84  // part
85 
86  v_ += 0.5*trackTime*a_;
87 
88  pi_ += 0.5*trackTime*tau_;
89  }
90  else if (td.part() == trackingData::tpLinearTrack)
91  {
92  // Leapfrog tracking part
93 
94  while (td.keepParticle && td.sendToProc == -1 && stepFraction() < 1)
95  {
96  const scalar f = 1 - stepFraction();
97  trackToAndHitFace(f*trackTime*v_, f, cloud, td);
98  }
99  }
100  else if (td.part() == trackingData::tpRotationalTrack)
101  {
102  // Leapfrog orientation adjustment, carried out before force calculation
103  // but after tracking stage, i.e. rotation carried once linear motion
104  // complete.
105 
106  if (!constProps.pointMolecule())
107  {
108  const diagTensor& momentOfInertia(constProps.momentOfInertia());
109 
110  tensor R;
111 
112  if (!constProps.linearMolecule())
113  {
114  R = rotationTensorX(0.5*trackTime*pi_.x()/momentOfInertia.xx());
115  pi_ = pi_ & R;
116  Q_ = Q_ & R;
117  }
118 
119  R = rotationTensorY(0.5*trackTime*pi_.y()/momentOfInertia.yy());
120  pi_ = pi_ & R;
121  Q_ = Q_ & R;
122 
123  R = rotationTensorZ(trackTime*pi_.z()/momentOfInertia.zz());
124  pi_ = pi_ & R;
125  Q_ = Q_ & R;
126 
127  R = rotationTensorY(0.5*trackTime*pi_.y()/momentOfInertia.yy());
128  pi_ = pi_ & R;
129  Q_ = Q_ & R;
130 
131  if (!constProps.linearMolecule())
132  {
133  R = rotationTensorX(0.5*trackTime*pi_.x()/momentOfInertia.xx());
134  pi_ = pi_ & R;
135  Q_ = Q_ & R;
136  }
137  }
138 
139  setSitePositions(td.mesh, constProps);
140  }
141  else if (td.part() == trackingData::tpVelocityHalfStep1)
142  {
143  // Second leapfrog velocity adjust part, required after tracking+force
144  // part
145 
146  scalar m = constProps.mass();
147 
148  a_ = Zero;
149 
150  tau_ = Zero;
151 
152  forAll(siteForces_, s)
153  {
154  const vector& f = siteForces_[s];
155 
156  a_ += f/m;
157 
158  tau_ += (constProps.siteReferencePositions()[s] ^ (Q_.T() & f));
159  }
160 
161  v_ += 0.5*trackTime*a_;
162 
163  pi_ += 0.5*trackTime*tau_;
164 
165  if (constProps.pointMolecule())
166  {
167  tau_ = Zero;
168 
169  pi_ = Zero;
170  }
171 
172  if (constProps.linearMolecule())
173  {
174  tau_.x() = 0.0;
175 
176  pi_.x() = 0.0;
177  }
178  }
179  else
180  {
182  << td.part() << " is an invalid part of the integration method."
183  << abort(FatalError);
184  }
185 
186  return td.keepParticle;
187 }
188 
189 
191 {
193 
194  Q_ = transform.T() & Q_;
195 
196  v_ = transform.transform(v_);
197 
198  a_ = transform.transform(a_);
199 
200  pi_ = Q_.T() & transform.transform(Q_ & pi_);
201 
202  tau_ = Q_.T() & transform.transform(Q_ & tau_);
203 
204  rf_ = transform.transform(rf_);
205 
206  transform.transformList(siteForces_);
207 
208  sitePositions_ = transform.transformPosition(vectorField(sitePositions_));
209 
210  if (special_ == SPECIAL_TETHERED)
211  {
212  specialPosition_ = transform.transformPosition(specialPosition_);
213  }
214 }
215 
216 
218 (
219  const polyMesh& mesh,
220  const constantProperties& constProps
221 )
222 {
223  sitePositions_ =
224  position(mesh)
225  + (Q_ & constProps.siteReferencePositions());
226 }
227 
228 
230 {
231  sitePositions_.setSize(size);
232 
233  siteForces_.setSize(size);
234 }
235 
236 
238 {
239  const vector nw = normal(td.mesh);
240  const scalar vn = v_ & nw;
241 
242  // Specular reflection
243  if (vn > 0)
244  {
245  v_ -= 2*vn*nw;
246  }
247 }
248 
249 
250 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
scalar deltaTValue() const
Return time step value.
Definition: TimeStateI.H:34
A cloud is a collection of lagrangian particles.
Definition: cloud.H:55
friend dimensionSet transform(const dimensionSet &)
Class to hold molecule constant properties.
Definition: molecule.H:91
const diagTensor & momentOfInertia() const
Definition: moleculeI.H:417
const Field< vector > & siteReferencePositions() const
Definition: moleculeI.H:332
Class used to pass tracking data to the trackToFace function.
Definition: molecule.H:169
trackPart part() const
Return the part of the tracking operation taking place.
Definition: molecule.H:203
void setSiteSizes(label size)
Definition: molecule.C:229
void setSitePositions(const polyMesh &mesh, const constantProperties &constProps)
Definition: molecule.C:218
bool move(moleculeCloud &, trackingData &)
Definition: molecule.C:69
virtual void transformProperties(const transformer &)
Transform the physical properties of the particle.
Definition: molecule.C:190
void hitWallPatch(moleculeCloud &cloud, trackingData &td)
Overridable function to handle the particle hitting a wallPatch.
Definition: molecule.C:237
Calculates the inertia tensor and principal axes and moments of a polyhedra/cells/triSurfaces....
const Time & time() const
Return time.
label sendToProc
Processor to send the particle to. -1 indicates that this.
Definition: particle.H:113
bool keepParticle
Flag to indicate whether to keep particle (false = delete)
Definition: particle.H:109
const polyMesh & mesh
Reference to the mesh.
Definition: particle.H:106
virtual void transformProperties(const transformer &)
Transform the physical properties of the particle.
Definition: particle.C:990
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:80
Vector-tensor class used to perform translations, rotations and scaling operations in 3D space.
Definition: transformer.H:84
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:334
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.name(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
static const zero Zero
Definition: zero.H:97
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
Tensor< scalar > tensor
Tensor of scalars.
Definition: tensor.H:51
errorManip< error > abort(error &err)
Definition: errorManip.H:131
dimensionedScalar sin(const dimensionedScalar &ds)
dimensionSet transform(const dimensionSet &)
Definition: dimensionSet.C:504
Field< vector > vectorField
Specialisation of Field<T> for vector.
static scalar R(const scalar a, const scalar x)
Definition: invIncGamma.C:102
error FatalError
dimensionedScalar cos(const dimensionedScalar &ds)
labelList f(nPoints)