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