src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 "KinematicParcel.H"
  28
  29 // * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
  30
  31 template<class ParcelType>
  32 template<class TrackData>
  33 void Foam::KinematicParcel<ParcelType>::setCellValues
  34 (
  35     TrackData& td,
  36     const scalar dt,
  37     const label cellI
  38 )
  39 {
  40     rhoc_ = td.rhoInterp().interpolate(this->position(), cellI);
  41     if (rhoc_ < td.constProps().rhoMin())
  42     {
  43         WarningIn
  44         (
  45             "void Foam::KinematicParcel<ParcelType>::setCellValues"
  46             "("
  47                 "TrackData&, "
  48                 "const scalar, "
  49                 "const label"
  50             ")"
  51         )   << "Limiting observed density in cell " << cellI << " to "
  52             << td.constProps().rhoMin() <<  nl << endl;
  53
  54         rhoc_ = td.constProps().rhoMin();
  55     }
  56
  57     Uc_ = td.UInterp().interpolate(this->position(), cellI);
  58
  59     muc_ = td.muInterp().interpolate(this->position(), cellI);
  60
  61     // Apply dispersion components to carrier phase velocity
  62     Uc_ = td.cloud().dispersion().update
  63     (
  64         dt,
  65         cellI,
  66         U_,
  67         Uc_,
  68         UTurb_,
  69         tTurb_
  70     );
  71 }
  72
  73
  74 template<class ParcelType>
  75 template<class TrackData>
  76 void Foam::KinematicParcel<ParcelType>::cellValueSourceCorrection
  77 (
  78     TrackData& td,
  79     const scalar dt,
  80     const label cellI
  81 )
  82 {
  83     Uc_ += td.cloud().UTrans()[cellI]/massCell(cellI);
  84 }
  85
  86
  87 template<class ParcelType>
  88 template<class TrackData>
  89 void Foam::KinematicParcel<ParcelType>::calc
  90 (
  91     TrackData& td,
  92     const scalar dt,
  93     const label cellI
  94 )
  95 {
  96     // Define local properties at beginning of time step
  97     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  98     const scalar np0 = nParticle_;
  99     const scalar d0 = d_;
 100     const vector U0 = U_;
 101     const scalar rho0 = rho_;
 102     const scalar mass0 = mass();
 103
 104     // Reynolds number
 105     const scalar Re = this->Re(U0, d0, rhoc_, muc_);
 106
 107
 108     // Sources
 109     //~~~~~~~~
 110
 111     // Explicit momentum source for particle
 112     vector Su = vector::zero;
 113
 114     // Momentum transfer from the particle to the carrier phase
 115     vector dUTrans = vector::zero;
 116
 117
 118     // Motion
 119     // ~~~~~~
 120
 121     // Calculate new particle velocity
 122     vector U1 =
 123         calcVelocity(td, dt, cellI, Re, muc_, d0, U0, rho0, mass0, Su, dUTrans);
 124
 125
 126     // Accumulate carrier phase source terms
 127     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 128     if (td.cloud().coupled())
 129     {
 130         // Update momentum transfer
 131         td.cloud().UTrans()[cellI] += np0*dUTrans;
 132     }
 133
 134
 135     // Set new particle properties
 136     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 137     U_ = U1;
 138 }
 139
 140
 141 template<class ParcelType>
 142 template<class TrackData>
 143 const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
 144 (
 145     TrackData& td,
 146     const scalar dt,
 147     const label cellI,
 148     const scalar Re,
 149     const scalar mu,
 150     const scalar d,
 151     const vector& U,
 152     const scalar rho,
 153     const scalar mass,
 154     const vector& Su,
 155     vector& dUTrans
 156 ) const
 157 {
 158     const polyMesh& mesh = this->cloud().pMesh();
 159
 160     // Momentum transfer coefficient
 161     const scalar utc = td.cloud().drag().utc(Re, d, mu) + ROOTVSMALL;
 162
 163     // Momentum source due to particle forces
 164     const vector FCoupled =
 165         mass*td.cloud().forces().calcCoupled(cellI, dt, rhoc_, rho, Uc_, U);
 166     const vector FNonCoupled =
 167         mass*td.cloud().forces().calcNonCoupled(cellI, dt, rhoc_, rho, Uc_, U);
 168
 169
 170     // New particle velocity
 171     //~~~~~~~~~~~~~~~~~~~~~~
 172
 173     // Update velocity - treat as 3-D
 174     const scalar As = this->areaS(d);
 175     const vector ap = Uc_ + (FCoupled + FNonCoupled + Su)/(utc*As);
 176     const scalar bp = 6.0*utc/(rho*d);
 177
 178     IntegrationScheme<vector>::integrationResult Ures =
 179         td.cloud().UIntegrator().integrate(U, dt, ap, bp);
 180
 181     vector Unew = Ures.value();
 182
 183     dUTrans += dt*(utc*As*(Ures.average() - Uc_) - FCoupled);
 184
 185     // Apply correction to velocity and dUTrans for reduced-D cases
 186     meshTools::constrainDirection(mesh, mesh.solutionD(), Unew);
 187     meshTools::constrainDirection(mesh, mesh.solutionD(), dUTrans);
 188
 189     return Unew;
 190 }
 191
 192
 193 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 194
 195 template <class ParcelType>
 196 Foam::KinematicParcel<ParcelType>::KinematicParcel
 197 (
 198     const KinematicParcel<ParcelType>& p
 199 )
 200 :
 201     Particle<ParcelType>(p),
 202     typeId_(p.typeId_),
 203     nParticle_(p.nParticle_),
 204     d_(p.d_),
 205     U_(p.U_),
 206     rho_(p.rho_),
 207     tTurb_(p.tTurb_),
 208     UTurb_(p.UTurb_),
 209     rhoc_(p.rhoc_),
 210     Uc_(p.Uc_),
 211     muc_(p.muc_)
 212 {}
 213
 214
 215 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 216
 217 template<class ParcelType>
 218 template<class TrackData>
 219 bool Foam::KinematicParcel<ParcelType>::move(TrackData& td)
 220 {
 221     ParcelType& p = static_cast<ParcelType&>(*this);
 222
 223     td.switchProcessor = false;
 224     td.keepParticle = true;
 225
 226     const polyMesh& mesh = td.cloud().pMesh();
 227     const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
 228
 229     const scalar deltaT = mesh.time().deltaTValue();
 230     scalar tEnd = (1.0 - p.stepFraction())*deltaT;
 231     const scalar dtMax = tEnd;
 232
 233     while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
 234     {
 235         // Apply correction to position for reduced-D cases
 236         meshTools::constrainToMeshCentre(mesh, p.position());
 237
 238         // Set the Lagrangian time-step
 239         scalar dt = min(dtMax, tEnd);
 240
 241         // Remember which cell the Parcel is in since this will change if a
 242         // face is hit
 243         label cellI = p.cell();
 244
 245         dt *= p.trackToFace(p.position() + dt*U_, td);
 246
 247         tEnd -= dt;
 248         p.stepFraction() = 1.0 - tEnd/deltaT;
 249
 250         // Avoid problems with extremely small timesteps
 251         if (dt > ROOTVSMALL)
 252         {
 253             // Update cell based properties
 254             p.setCellValues(td, dt, cellI);
 255
 256             if (td.cloud().cellValueSourceCorrection())
 257             {
 258                 p.cellValueSourceCorrection(td, dt, cellI);
 259             }
 260
 261             p.calc(td, dt, cellI);
 262         }
 263
 264         if (p.onBoundary() && td.keepParticle)
 265         {
 266             if (isA<processorPolyPatch>(pbMesh[p.patch(p.face())]))
 267             {
 268                 td.switchProcessor = true;
 269             }
 270         }
 271     }
 272
 273     return td.keepParticle;
 274 }
 275
 276
 277 template<class ParcelType>
 278 template<class TrackData>
 279 bool Foam::KinematicParcel<ParcelType>::hitPatch
 280 (
 281     const polyPatch& pp,
 282     TrackData& td,
 283     const label patchI
 284 )
 285 {
 286     ParcelType& p = static_cast<ParcelType&>(*this);
 287     td.cloud().postProcessing().postPatch(p, patchI);
 288
 289     return td.cloud().patchInteraction().correct
 290     (
 291         pp,
 292         this->face(),
 293         td.keepParticle,
 294         U_
 295     );
 296 }
 297
 298
 299 template<class ParcelType>
 300 bool Foam::KinematicParcel<ParcelType>::hitPatch
 301 (
 302     const polyPatch& pp,
 303     int& td,
 304     const label patchI
 305 )
 306 {
 307     return false;
 308 }
 309
 310
 311 template<class ParcelType>
 312 template<class TrackData>
 313 void Foam::KinematicParcel<ParcelType>::hitProcessorPatch
 314 (
 315     const processorPolyPatch&,
 316     TrackData& td
 317 )
 318 {
 319     td.switchProcessor = true;
 320 }
 321
 322
 323 template<class ParcelType>
 324 void Foam::KinematicParcel<ParcelType>::hitProcessorPatch
 325 (
 326     const processorPolyPatch&,
 327     int&
 328 )
 329 {}
 330
 331
 332 template<class ParcelType>
 333 template<class TrackData>
 334 void Foam::KinematicParcel<ParcelType>::hitWallPatch
 335 (
 336     const wallPolyPatch& wpp,
 337     TrackData& td
 338 )
 339 {
 340     // Wall interactions handled by generic hitPatch function
 341 }
 342
 343
 344 template<class ParcelType>
 345 void Foam::KinematicParcel<ParcelType>::hitWallPatch
 346 (
 347     const wallPolyPatch&,
 348     int&
 349 )
 350 {}
 351
 352
 353 template<class ParcelType>
 354 template<class TrackData>
 355 void Foam::KinematicParcel<ParcelType>::hitPatch
 356 (
 357     const polyPatch&,
 358     TrackData& td
 359 )
 360 {
 361     td.keepParticle = false;
 362 }
 363
 364
 365 template<class ParcelType>
 366 void Foam::KinematicParcel<ParcelType>::hitPatch(const polyPatch&, int&)
 367 {}
 368
 369
 370 template<class ParcelType>
 371 void Foam::KinematicParcel<ParcelType>::transformProperties(const tensor& T)
 372 {
 373     Particle<ParcelType>::transformProperties(T);
 374     U_ = transform(T, U_);
 375 }
 376
 377
 378 template<class ParcelType>
 379 void Foam::KinematicParcel<ParcelType>::transformProperties
 380 (
 381     const vector& separation
 382 )
 383 {
 384     Particle<ParcelType>::transformProperties(separation);
 385 }
 386
 387
 388 // * * * * * * * * * * * * * * IOStream operators  * * * * * * * * * * * * * //
 389
 390 #include "KinematicParcelIO.C"
 391
 392 // ************************************************************************* //
 393