DSMCParcel.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 "DSMCParcel.H"
27 #include "meshTools.H"
28 
29 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
30 
31 template<class ParcelType>
32 template<class TrackCloudType>
34 (
35  TrackCloudType& cloud,
36  trackingData& td,
37  const scalar trackTime
38 )
39 {
40  typename TrackCloudType::parcelType& p =
41  static_cast<typename TrackCloudType::parcelType&>(*this);
42 
43  td.switchProcessor = false;
44  td.keepParticle = true;
45 
46  const polyMesh& mesh = cloud.pMesh();
47 
48  // For reduced-D cases, the velocity used to track needs to be
49  // constrained, but the actual U_ of the parcel must not be
50  // altered or used, as it is altered by patch interactions an
51  // needs to retain its 3D value for collision purposes.
52  vector Utracking = U_;
53 
54  while (td.keepParticle && !td.switchProcessor && p.stepFraction() < 1)
55  {
56  Utracking = U_;
57 
58  // Apply correction to velocity to constrain tracking for
59  // reduced-D cases
60  meshTools::constrainDirection(mesh, mesh.solutionD(), Utracking);
61 
62  // Deviation from the mesh centre for reduced-D cases
63  const vector d = p.deviationFromMeshCentre();
64 
65  const scalar f = 1 - p.stepFraction();
66  p.trackToAndHitFace(f*trackTime*Utracking - d, f, cloud, td);
67  }
68 
69  return td.keepParticle;
70 }
71 
72 
73 template<class ParcelType>
74 template<class TrackCloudType>
76 {
77  return false;
78 }
79 
80 
81 template<class ParcelType>
82 template<class TrackCloudType>
84 (
85  TrackCloudType&,
86  trackingData& td
87 )
88 {
89  td.switchProcessor = true;
90 }
91 
92 
93 template<class ParcelType>
94 template<class TrackCloudType>
96 (
97  TrackCloudType& cloud,
99 )
100 {
101  const label wppIndex = this->patch();
102 
103  const wallPolyPatch& wpp =
104  static_cast<const wallPolyPatch&>
105  (
106  this->mesh().boundaryMesh()[wppIndex]
107  );
108 
109  const label wppLocalFace = wpp.whichFace(this->face());
110 
111  const scalar fA = wpp.magFaceAreas()[wppLocalFace];
112 
113  const scalar deltaT = cloud.pMesh().time().deltaTValue();
114 
115  const constantProperties& constProps(cloud.constProps(typeId_));
116 
117  scalar m = constProps.mass();
118 
119  vector nw = wpp.faceAreas()[wppLocalFace];
120  nw /= mag(nw);
121 
122  scalar U_dot_nw = U_ & nw;
123 
124  vector Ut = U_ - U_dot_nw*nw;
125 
126  scalar invMagUnfA = 1/max(mag(U_dot_nw)*fA, vSmall);
127 
128  cloud.rhoNBF()[wppIndex][wppLocalFace] += invMagUnfA;
129 
130  cloud.rhoMBF()[wppIndex][wppLocalFace] += m*invMagUnfA;
131 
132  cloud.linearKEBF()[wppIndex][wppLocalFace] +=
133  0.5*m*(U_ & U_)*invMagUnfA;
134 
135  cloud.internalEBF()[wppIndex][wppLocalFace] += Ei_*invMagUnfA;
136 
137  cloud.iDofBF()[wppIndex][wppLocalFace] +=
138  constProps.internalDegreesOfFreedom()*invMagUnfA;
139 
140  cloud.momentumBF()[wppIndex][wppLocalFace] += m*Ut*invMagUnfA;
141 
142  // pre-interaction energy
143  scalar preIE = 0.5*m*(U_ & U_) + Ei_;
144 
145  // pre-interaction momentum
146  vector preIMom = m*U_;
147 
148  cloud.wallInteraction().correct(*this);
149 
150  U_dot_nw = U_ & nw;
151 
152  Ut = U_ - U_dot_nw*nw;
153 
154  invMagUnfA = 1/max(mag(U_dot_nw)*fA, vSmall);
155 
156  cloud.rhoNBF()[wppIndex][wppLocalFace] += invMagUnfA;
157 
158  cloud.rhoMBF()[wppIndex][wppLocalFace] += m*invMagUnfA;
159 
160  cloud.linearKEBF()[wppIndex][wppLocalFace] +=
161  0.5*m*(U_ & U_)*invMagUnfA;
162 
163  cloud.internalEBF()[wppIndex][wppLocalFace] += Ei_*invMagUnfA;
164 
165  cloud.iDofBF()[wppIndex][wppLocalFace] +=
166  constProps.internalDegreesOfFreedom()*invMagUnfA;
167 
168  cloud.momentumBF()[wppIndex][wppLocalFace] += m*Ut*invMagUnfA;
169 
170  // post-interaction energy
171  scalar postIE = 0.5*m*(U_ & U_) + Ei_;
172 
173  // post-interaction momentum
174  vector postIMom = m*U_;
175 
176  scalar deltaQ = cloud.nParticle()*(preIE - postIE)/(deltaT*fA);
177 
178  vector deltaFD = cloud.nParticle()*(preIMom - postIMom)/(deltaT*fA);
179 
180  cloud.qBF()[wppIndex][wppLocalFace] += deltaQ;
181 
182  cloud.fDBF()[wppIndex][wppLocalFace] += deltaFD;
183 }
184 
185 
186 template<class ParcelType>
188 (
189  const transformer& transform
190 )
191 {
192  ParcelType::transformProperties(transform);
193  U_ = transform.transform(U_);
194 }
195 
196 
197 // * * * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * //
198 
199 #include "DSMCParcelIO.C"
200 
201 // ************************************************************************* //
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
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
const Vector< label > & solutionD() const
Return the vector of solved-for directions in mesh.
Definition: polyMesh.C:831
void hitWallPatch(TrackCloudType &, trackingData &)
Overridable function to handle the particle hitting a wallPatch.
Definition: DSMCParcel.C:96
dynamicFvMesh & mesh
Class to hold DSMC particle constant properties.
Definition: DSMCParcel.H:80
Foam::wallPolyPatch.
Definition: wallPolyPatch.H:48
scalar mass() const
Return const access to the particle mass [kg].
Definition: DSMCParcelI.H:96
const vectorField::subField faceAreas() const
Return face areas.
Definition: polyPatch.C:303
labelList f(nPoints)
void hitProcessorPatch(TrackCloudType &, trackingData &)
Overridable function to handle the particle hitting a.
Definition: DSMCParcel.C:84
bool move(TrackCloudType &cloud, trackingData &td, const scalar trackTime)
Move the parcel.
Definition: DSMCParcel.C:34
const scalarField::subField magFaceAreas() const
Return face area magnitudes.
Definition: polyPatch.C:309
virtual void transformProperties(const transformer &)
Transform the physical properties of the particle.
Definition: DSMCParcel.C:188
bool hitPatch(TrackCloudType &, trackingData &)
Overridable function to handle the particle hitting a patch.
Definition: DSMCParcel.C:75
ParcelType::trackingData trackingData
Use base tracking data.
Definition: DSMCParcel.H:130
dimensioned< scalar > mag(const dimensioned< Type > &)
Type transform(const Type &) const
Transform the given type.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
void constrainDirection(const polyMesh &mesh, const Vector< label > &dirs, vector &d)
Set the constrained components of directions/velocity to zero.
Definition: meshTools.C:671
volScalarField & p
label nw
Definition: createFields.H:12
label whichFace(const label l) const
Return label of face in patch from global face label.
Definition: polyPatch.H:389