FreeFOAM The Cross-Platform CFD Toolkit
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) 2008-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 \*----------------------------------------------------------------------------*/
25 
26 #include <molecule/moleculeCloud.H>
27 #include "molecule.H"
28 #include <OpenFOAM/Random.H>
29 #include <OpenFOAM/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 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
67 
69 (
70  moleculeCloud& molCloud,
71  label part
72 )
73 :
75  molCloud_(molCloud),
76  part_(part)
77 {}
78 
79 
80 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
81 
82 
84 {
85  td.switchProcessor = false;
86  td.keepParticle = true;
87 
88  const constantProperties& constProps(td.molCloud().constProps(id_));
89 
90  scalar deltaT = cloud().pMesh().time().deltaT().value();
91 
92  if (td.part() == 0)
93  {
94  // First leapfrog velocity adjust part, required before tracking+force
95  // part
96 
97  v_ += 0.5*deltaT*a_;
98 
99  pi_ += 0.5*deltaT*tau_;
100  }
101  else if (td.part() == 1)
102  {
103  // Leapfrog tracking part
104 
105  scalar tEnd = (1.0 - stepFraction())*deltaT;
106  scalar dtMax = tEnd;
107 
108  while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
109  {
110  // set the lagrangian time-step
111  scalar dt = min(dtMax, tEnd);
112 
113  dt *= trackToFace(position() + dt*v_, td);
114 
115  tEnd -= dt;
116  stepFraction() = 1.0 - tEnd/deltaT;
117  }
118  }
119  else if (td.part() == 2)
120  {
121  // Leapfrog orientation adjustment, carried out before force calculation
122  // but after tracking stage, i.e. rotation carried once linear motion
123  // complete.
124 
125  if (!constProps.pointMolecule())
126  {
127  const diagTensor& momentOfInertia(constProps.momentOfInertia());
128 
129  tensor R;
130 
131  if (!constProps.linearMolecule())
132  {
133  R = rotationTensorX(0.5*deltaT*pi_.x()/momentOfInertia.xx());
134  pi_ = pi_ & R;
135  Q_ = Q_ & R;
136  }
137 
138  R = rotationTensorY(0.5*deltaT*pi_.y()/momentOfInertia.yy());
139  pi_ = pi_ & R;
140  Q_ = Q_ & R;
141 
142  R = rotationTensorZ(deltaT*pi_.z()/momentOfInertia.zz());
143  pi_ = pi_ & R;
144  Q_ = Q_ & R;
145 
146  R = rotationTensorY(0.5*deltaT*pi_.y()/momentOfInertia.yy());
147  pi_ = pi_ & R;
148  Q_ = Q_ & R;
149 
150  if (!constProps.linearMolecule())
151  {
152  R = rotationTensorX(0.5*deltaT*pi_.x()/momentOfInertia.xx());
153  pi_ = pi_ & R;
154  Q_ = Q_ & R;
155  }
156  }
157 
158  setSitePositions(constProps);
159  }
160  else if (td.part() == 3)
161  {
162  // Second leapfrog velocity adjust part, required after tracking+force
163  // part
164 
165  scalar m = constProps.mass();
166 
167  a_ = vector::zero;
168 
169  tau_ = vector::zero;
170 
171  forAll(siteForces_, s)
172  {
173  const vector& f = siteForces_[s];
174 
175  a_ += f/m;
176 
177  tau_ += (constProps.siteReferencePositions()[s] ^ (Q_.T() & f));
178  }
179 
180  v_ += 0.5*deltaT*a_;
181 
182  pi_ += 0.5*deltaT*tau_;
183 
184  if (constProps.pointMolecule())
185  {
186  tau_ = vector::zero;
187 
188  pi_ = vector::zero;
189  }
190 
191  if (constProps.linearMolecule())
192  {
193  tau_.x() = 0.0;
194 
195  pi_.x() = 0.0;
196  }
197  }
198  else
199  {
200  FatalErrorIn("molecule::move(molecule::trackData& td)") << nl
201  << td.part()
202  << " is an invalid part of the integration method."
203  << abort(FatalError);
204  }
205 
206  return td.keepParticle;
207 }
208 
209 
211 {
212  Q_ = T & Q_;
213 
214  sitePositions_ = position_ + (T & (sitePositions_ - position_));
215 
216  siteForces_ = T & siteForces_;
217 }
218 
219 
221 {
222  if (special_ == SPECIAL_TETHERED)
223  {
224  specialPosition_ += separation;
225  }
226 }
227 
228 
230 {
231  sitePositions_ = position_ + (Q_ & constProps.siteReferencePositions());
232 }
233 
234 
236 {
237  sitePositions_.setSize(size);
238 
239  siteForces_.setSize(size);
240 }
241 
242 
244 (
245  const polyPatch&,
247  const label
248 )
249 {
250  return false;
251 }
252 
253 
255 (
256  const polyPatch&,
257  int&,
258  const label
259 )
260 {
261  return false;
262 }
263 
264 
266 (
267  const processorPolyPatch&,
269 )
270 {
271  td.switchProcessor = true;
272 }
273 
274 
276 (
277  const processorPolyPatch&,
278  int&
279 )
280 {}
281 
282 
284 (
285  const wallPolyPatch& wpp,
287 )
288 {
289  vector nw = wpp.faceAreas()[wpp.whichFace(face())];
290  nw /= mag(nw);
291 
292  scalar vn = v_ & nw;
293 
294  // Specular reflection
295  if (vn > 0)
296  {
297  v_ -= 2*vn*nw;
298  }
299 }
300 
301 
303 (
304  const wallPolyPatch&,
305  int&
306 )
307 {}
308 
309 
311 (
312  const polyPatch&,
314 )
315 {
316  td.keepParticle = false;
317 }
318 
319 
321 (
322  const polyPatch&,
323  int&
324 )
325 {}
326 
327 
328 // ************************ vim: set sw=4 sts=4 et: ************************ //