src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C

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