FreeFOAM The Cross-Platform CFD Toolkit
DsmcParcel_.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) 2009-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 "DsmcParcel_.H"
27 #include <meshTools/meshTools.H>
28 
29 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
30 
31 template<class ParcelType>
32 template<class TrackData>
34 (
35  TrackData& td
36 )
37 {
38  ParcelType& p = static_cast<ParcelType&>(*this);
39 
40  td.switchProcessor = false;
41  td.keepParticle = true;
42 
43  const polyMesh& mesh = td.cloud().pMesh();
44  const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
45 
46  const scalar deltaT = mesh.time().deltaTValue();
47  scalar tEnd = (1.0 - p.stepFraction())*deltaT;
48  const scalar dtMax = tEnd;
49 
50  // For reduced-D cases, the velocity used to track needs to be
51  // constrained, but the actual U_ of the parcel must not be
52  // altered or used, as it is altered by patch interactions an
53  // needs to retain its 3D value for collision purposes.
54  vector Utracking = U_;
55 
56  while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
57  {
58  // Apply correction to position for reduced-D cases
59  meshTools::constrainToMeshCentre(mesh, p.position());
60 
61  Utracking = U_;
62 
63  // Apply correction to velocity to constrain tracking for
64  // reduced-D cases
65  meshTools::constrainDirection(mesh, mesh.solutionD(), Utracking);
66 
67  // Set the Lagrangian time-step
68  scalar dt = min(dtMax, tEnd);
69 
70  dt *= p.trackToFace(p.position() + dt*Utracking, td);
71 
72  tEnd -= dt;
73 
74  p.stepFraction() = 1.0 - tEnd/deltaT;
75 
76  if (p.onBoundary() && td.keepParticle)
77  {
78  if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())]))
79  {
80  td.switchProcessor = true;
81  }
82  }
83  }
84 
85  return td.keepParticle;
86 }
87 
88 
89 template<class ParcelType>
90 template<class TrackData>
92 (
93  const polyPatch&,
94  TrackData& td,
95  const label patchI
96 )
97 {
98  return false;
99 }
100 
101 
102 template<class ParcelType>
103 template<class TrackData>
105 (
106  const processorPolyPatch&,
107  TrackData& td
108 )
109 {
110  td.switchProcessor = true;
111 }
112 
113 
114 template<class ParcelType>
116 (
117  const processorPolyPatch&,
118  int&
119 )
120 {}
121 
122 
123 template<class ParcelType>
124 template<class TrackData>
126 (
127  const wallPolyPatch& wpp,
128  TrackData& td
129 )
130 {
131  label wppIndex = wpp.index();
132 
133  label wppLocalFace = wpp.whichFace(this->face());
134 
135  const scalar fA = mag(wpp.faceAreas()[wppLocalFace]);
136 
137  const scalar deltaT = td.cloud().pMesh().time().deltaTValue();
138 
139  const constantProperties& constProps(td.cloud().constProps(typeId_));
140 
141  scalar m = constProps.mass();
142 
143  vector nw = wpp.faceAreas()[wppLocalFace];
144  nw /= mag(nw);
145 
146  scalar U_dot_nw = U_ & nw;
147 
148  vector Ut = U_ - U_dot_nw*nw;
149 
150  scalar invMagUnfA = 1/max(mag(U_dot_nw)*fA, VSMALL);
151 
152  td.cloud().rhoNBF()[wppIndex][wppLocalFace] += invMagUnfA;
153 
154  td.cloud().rhoMBF()[wppIndex][wppLocalFace] += m*invMagUnfA;
155 
156  td.cloud().linearKEBF()[wppIndex][wppLocalFace] +=
157  0.5*m*(U_ & U_)*invMagUnfA;
158 
159  td.cloud().internalEBF()[wppIndex][wppLocalFace] += Ei_*invMagUnfA;
160 
161  td.cloud().iDofBF()[wppIndex][wppLocalFace] +=
162  constProps.internalDegreesOfFreedom()*invMagUnfA;
163 
164  td.cloud().momentumBF()[wppIndex][wppLocalFace] += m*Ut*invMagUnfA;
165 
166  // pre-interaction energy
167  scalar preIE = 0.5*m*(U_ & U_) + Ei_;
168 
169  // pre-interaction momentum
170  vector preIMom = m*U_;
171 
172  td.cloud().wallInteraction().correct
173  (
174  wpp,
175  this->face(),
176  U_,
177  Ei_,
178  typeId_
179  );
180 
181  U_dot_nw = U_ & nw;
182 
183  Ut = U_ - U_dot_nw*nw;
184 
185  invMagUnfA = 1/max(mag(U_dot_nw)*fA, VSMALL);
186 
187  td.cloud().rhoNBF()[wppIndex][wppLocalFace] += invMagUnfA;
188 
189  td.cloud().rhoMBF()[wppIndex][wppLocalFace] += m*invMagUnfA;
190 
191  td.cloud().linearKEBF()[wppIndex][wppLocalFace] +=
192  0.5*m*(U_ & U_)*invMagUnfA;
193 
194  td.cloud().internalEBF()[wppIndex][wppLocalFace] += Ei_*invMagUnfA;
195 
196  td.cloud().iDofBF()[wppIndex][wppLocalFace] +=
197  constProps.internalDegreesOfFreedom()*invMagUnfA;
198 
199  td.cloud().momentumBF()[wppIndex][wppLocalFace] += m*Ut*invMagUnfA;
200 
201  // post-interaction energy
202  scalar postIE = 0.5*m*(U_ & U_) + Ei_;
203 
204  // post-interaction momentum
205  vector postIMom = m*U_;
206 
207  scalar deltaQ = td.cloud().nParticle()*(preIE - postIE)/(deltaT*fA);
208 
209  vector deltaFD = td.cloud().nParticle()*(preIMom - postIMom)/(deltaT*fA);
210 
211  td.cloud().qBF()[wppIndex][wppLocalFace] += deltaQ;
212 
213  td.cloud().fDBF()[wppIndex][wppLocalFace] += deltaFD;
214 
215 }
216 
217 
218 template<class ParcelType>
220 (
221  const wallPolyPatch&,
222  int&
223 )
224 {}
225 
226 
227 template<class ParcelType>
228 template<class TrackData>
230 (
231  const polyPatch&,
232  TrackData& td
233 )
234 {
235  td.keepParticle = false;
236 }
237 
238 
239 template<class ParcelType>
241 (
242  const polyPatch&,
243  int&
244 )
245 {}
246 
247 
248 template<class ParcelType>
250 (
251  const tensor& T
252 )
253 {
255  U_ = transform(T, U_);
256 }
257 
258 
259 template<class ParcelType>
261 (
262  const vector& separation
263 )
264 {
266 }
267 
268 
269 // * * * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * //
270 
271 #include "DsmcParcelIO.C"
272 
273 // ************************ vim: set sw=4 sts=4 et: ************************ //