src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2008-2011 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 "fvmSup.H"
  27
  28 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  29
  30 template<class CloudType>
  31 inline const Foam::KinematicCloud<CloudType>&
  32 Foam::KinematicCloud<CloudType>::cloudCopy() const
  33 {
  34     return cloudCopyPtr_();
  35 }
  36
  37
  38 template<class CloudType>
  39 inline const Foam::fvMesh& Foam::KinematicCloud<CloudType>::mesh() const
  40 {
  41     return mesh_;
  42 }
  43
  44
  45 template<class CloudType>
  46 inline const Foam::IOdictionary&
  47 Foam::KinematicCloud<CloudType>::particleProperties() const
  48 {
  49     return particleProperties_;
  50 }
  51
  52
  53 template<class CloudType>
  54 inline const Foam::cloudSolution&
  55 Foam::KinematicCloud<CloudType>::solution() const
  56 {
  57     return solution_;
  58 }
  59
  60
  61 template<class CloudType>
  62 inline Foam::cloudSolution& Foam::KinematicCloud<CloudType>::solution()
  63 {
  64     return solution_;
  65 }
  66
  67
  68 template<class CloudType>
  69 inline const typename CloudType::particleType::constantProperties&
  70 Foam::KinematicCloud<CloudType>::constProps() const
  71 {
  72     return constProps_;
  73 }
  74
  75
  76 template<class CloudType>
  77 inline const Foam::dictionary&
  78 Foam::KinematicCloud<CloudType>::subModelProperties() const
  79 {
  80     return subModelProperties_;
  81 }
  82
  83
  84 template<class CloudType>
  85 inline const Foam::volScalarField& Foam::KinematicCloud<CloudType>::rho() const
  86 {
  87     return rho_;
  88 }
  89
  90
  91 template<class CloudType>
  92 inline const Foam::volVectorField& Foam::KinematicCloud<CloudType>::U() const
  93 {
  94     return U_;
  95 }
  96
  97
  98 template<class CloudType>
  99 inline const Foam::volScalarField& Foam::KinematicCloud<CloudType>::mu() const
 100 {
 101     return mu_;
 102 }
 103
 104
 105 template<class CloudType>
 106 inline const Foam::dimensionedVector& Foam::KinematicCloud<CloudType>::g() const
 107 {
 108     return g_;
 109 }
 110
 111
 112 template<class CloudType>
 113 inline Foam::scalar Foam::KinematicCloud<CloudType>::pAmbient() const
 114 {
 115     return pAmbient_;
 116 }
 117
 118
 119 template<class CloudType>
 120 inline Foam::scalar& Foam::KinematicCloud<CloudType>::pAmbient()
 121 {
 122     return pAmbient_;
 123 }
 124
 125
 126 template<class CloudType>
 127 //inline const typename CloudType::parcelType::forceType&
 128 inline const typename Foam::KinematicCloud<CloudType>::forceType&
 129 Foam::KinematicCloud<CloudType>::forces() const
 130 {
 131     return forces_;
 132 }
 133
 134
 135 template<class CloudType>
 136 inline typename Foam::KinematicCloud<CloudType>::functionType&
 137 Foam::KinematicCloud<CloudType>::functions()
 138 {
 139     return functions_;
 140 }
 141
 142
 143 template<class CloudType>
 144 inline const Foam::DispersionModel<Foam::KinematicCloud<CloudType> >&
 145 Foam::KinematicCloud<CloudType>::dispersion() const
 146 {
 147     return dispersionModel_;
 148 }
 149
 150
 151 template<class CloudType>
 152 inline Foam::DispersionModel<Foam::KinematicCloud<CloudType> >&
 153 Foam::KinematicCloud<CloudType>::dispersion()
 154 {
 155     return dispersionModel_();
 156 }
 157
 158
 159 template<class CloudType>
 160 inline const Foam::InjectionModel<Foam::KinematicCloud<CloudType> >&
 161 Foam::KinematicCloud<CloudType>::injection() const
 162 {
 163     return injectionModel_;
 164 }
 165
 166
 167 template<class CloudType>
 168 inline const Foam::PatchInteractionModel<Foam::KinematicCloud<CloudType> >&
 169 Foam::KinematicCloud<CloudType>::patchInteraction() const
 170 {
 171     return patchInteractionModel_;
 172 }
 173
 174
 175 template<class CloudType>
 176 inline Foam::PatchInteractionModel<Foam::KinematicCloud<CloudType> >&
 177 Foam::KinematicCloud<CloudType>::patchInteraction()
 178 {
 179     return patchInteractionModel_();
 180 }
 181
 182
 183 template<class CloudType>
 184 inline Foam::InjectionModel<Foam::KinematicCloud<CloudType> >&
 185 Foam::KinematicCloud<CloudType>::injection()
 186 {
 187     return injectionModel_();
 188 }
 189
 190
 191 template<class CloudType>
 192 inline const Foam::SurfaceFilmModel<Foam::KinematicCloud<CloudType> >&
 193 Foam::KinematicCloud<CloudType>::surfaceFilm() const
 194 {
 195     return surfaceFilmModel_();
 196 }
 197
 198
 199 template<class CloudType>
 200 inline Foam::SurfaceFilmModel<Foam::KinematicCloud<CloudType> >&
 201 Foam::KinematicCloud<CloudType>::surfaceFilm()
 202 {
 203     return surfaceFilmModel_();
 204 }
 205
 206
 207 template<class CloudType>
 208 inline const Foam::vectorIntegrationScheme&
 209 Foam::KinematicCloud<CloudType>::UIntegrator() const
 210 {
 211     return UIntegrator_;
 212 }
 213
 214
 215 template<class CloudType>
 216 inline Foam::scalar Foam::KinematicCloud<CloudType>::massInSystem() const
 217 {
 218     scalar sysMass = 0.0;
 219     forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
 220     {
 221          const parcelType& p = iter();
 222          sysMass += p.mass()*p.nParticle();
 223     }
 224
 225     return sysMass;
 226 }
 227
 228
 229 template<class CloudType>
 230 inline Foam::vector
 231 Foam::KinematicCloud<CloudType>::linearMomentumOfSystem() const
 232 {
 233     vector linearMomentum(vector::zero);
 234
 235     forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
 236     {
 237         const parcelType& p = iter();
 238
 239         linearMomentum += p.mass()*p.U();
 240     }
 241
 242     return linearMomentum;
 243 }
 244
 245
 246 template<class CloudType>
 247 inline Foam::scalar
 248 Foam::KinematicCloud<CloudType>::linearKineticEnergyOfSystem() const
 249 {
 250     scalar linearKineticEnergy = 0.0;
 251
 252     forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
 253     {
 254         const parcelType& p = iter();
 255
 256         linearKineticEnergy += 0.5*p.mass()*(p.U() & p.U());
 257     }
 258
 259     return linearKineticEnergy;
 260 }
 261
 262
 263 template<class CloudType>
 264 inline Foam::scalar
 265 Foam::KinematicCloud<CloudType>::rotationalKineticEnergyOfSystem() const
 266 {
 267     scalar rotationalKineticEnergy = 0.0;
 268
 269     forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
 270     {
 271         const parcelType& p = iter();
 272
 273         rotationalKineticEnergy +=
 274             0.5*p.momentOfInertia()*(p.omega() & p.omega());
 275     }
 276
 277     return rotationalKineticEnergy;
 278 }
 279
 280
 281 template<class CloudType>
 282 inline Foam::scalar Foam::KinematicCloud<CloudType>::Dij
 283 (
 284     const label i,
 285     const label j
 286 ) const
 287 {
 288     scalar si = 0.0;
 289     scalar sj = 0.0;
 290     forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
 291     {
 292         const parcelType& p = iter();
 293         si += p.nParticle()*pow(p.d(), i);
 294         sj += p.nParticle()*pow(p.d(), j);
 295     }
 296
 297     reduce(si, sumOp<scalar>());
 298     reduce(sj, sumOp<scalar>());
 299     sj = max(sj, VSMALL);
 300
 301     return si/sj;
 302 }
 303
 304
 305 template<class CloudType>
 306 inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration
 307 (
 308     const scalar& prc
 309 ) const
 310 {
 311     scalar distance = 0.0;
 312     scalar mTot = 0.0;
 313
 314     label np = this->size();
 315
 316     // arrays containing the parcels mass and
 317     // distance from injector in ascending order
 318     scalarField mass(np);
 319     scalarField dist(np);
 320
 321     if (np > 0)
 322     {
 323         label n = 0;
 324
 325         // first arrange the parcels in ascending order
 326         // the first parcel is closest to its injection position
 327         // and the last one is most far away.
 328         forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
 329         {
 330             const parcelType& p = iter();
 331             scalar mi = p.nParticle()*p.mass();
 332             scalar di = mag(p.position() - p.position0());
 333             mTot += mi;
 334
 335             // insert at the last place
 336             mass[n] = mi;
 337             dist[n] = di;
 338
 339             label i = 0;
 340             bool found = false;
 341
 342             // insert the parcel in the correct place
 343             // and move the others
 344             while ((i < n) && (!found))
 345             {
 346                 if (di < dist[i])
 347                 {
 348                     found = true;
 349                     for (label j=n; j>i; j--)
 350                     {
 351                         mass[j] = mass[j-1];
 352                         dist[j] = dist[j-1];
 353                     }
 354                     mass[i] = mi;
 355                     dist[i] = di;
 356                 }
 357                 i++;
 358             }
 359             n++;
 360         }
 361     }
 362
 363     reduce(mTot, sumOp<scalar>());
 364
 365     if (np > 0)
 366     {
 367         scalar mLimit = prc*mTot;
 368         scalar mOff = (1.0 - prc)*mTot;
 369
 370         if (np > 1)
 371         {
 372             // 'prc' is large enough that the parcel most far
 373             // away will be used, no need to loop...
 374             if (mLimit > mTot - mass[np-1])
 375             {
 376                 distance = dist[np-1];
 377             }
 378             else
 379             {
 380                 scalar mOffSum = 0.0;
 381                 label i = np;
 382
 383                 while ((mOffSum < mOff) && (i>0))
 384                 {
 385                     i--;
 386                     mOffSum += mass[i];
 387                 }
 388                 distance =
 389                     dist[i+1]
 390                   + (dist[i] - dist[i+1])*(mOffSum - mOff)
 391                    /mass[i+1] ;
 392             }
 393         }
 394         else
 395         {
 396             distance = dist[0];
 397         }
 398     }
 399
 400     reduce(distance, maxOp<scalar>());
 401
 402     return distance;
 403 }
 404
 405
 406 template<class CloudType>
 407 inline Foam::cachedRandom& Foam::KinematicCloud<CloudType>::rndGen()
 408 {
 409     return rndGen_;
 410 }
 411
 412
 413 template<class CloudType>
 414 inline Foam::List<Foam::DynamicList<typename CloudType::particleType*> >&
 415 Foam::KinematicCloud<CloudType>::cellOccupancy()
 416 {
 417     if (cellOccupancyPtr_.empty())
 418     {
 419         buildCellOccupancy();
 420     }
 421
 422     return cellOccupancyPtr_();
 423 }
 424
 425
 426 template<class CloudType>
 427 inline Foam::DimensionedField<Foam::vector, Foam::volMesh>&
 428 Foam::KinematicCloud<CloudType>::UTrans()
 429 {
 430     return UTrans_();
 431 }
 432
 433
 434 template<class CloudType>
 435 inline const Foam::DimensionedField<Foam::vector, Foam::volMesh>&
 436 Foam::KinematicCloud<CloudType>::UTrans() const
 437 {
 438     return UTrans_();
 439 }
 440
 441
 442 template<class CloudType>
 443 inline Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
 444 Foam::KinematicCloud<CloudType>::UCoeff()
 445 {
 446     return UCoeff_();
 447 }
 448
 449
 450 template<class CloudType>
 451 inline const Foam::DimensionedField<Foam::scalar, Foam::volMesh>&
 452 Foam::KinematicCloud<CloudType>::UCoeff() const
 453 {
 454     return UCoeff_();
 455 }
 456
 457
 458 template<class CloudType>
 459 inline Foam::tmp<Foam::fvVectorMatrix>
 460 Foam::KinematicCloud<CloudType>::SU(volVectorField& U) const
 461 {
 462     if (debug)
 463     {
 464         Info<< "UTrans min/max = " << min(UTrans()).value() << ", "
 465             << max(UTrans()).value() << nl
 466             << "UCoeff min/max = " << min(UCoeff()).value() << ", "
 467             << max(UCoeff()).value() << endl;
 468     }
 469
 470     if (solution_.coupled())
 471     {
 472         if (solution_.semiImplicit("U"))
 473         {
 474             const DimensionedField<scalar, volMesh>
 475                 Vdt(mesh_.V()*this->db().time().deltaT());
 476
 477             return UTrans()/Vdt - fvm::Sp(UCoeff()/Vdt, U) + UCoeff()/Vdt*U;
 478         }
 479         else
 480         {
 481             tmp<fvVectorMatrix> tfvm(new fvVectorMatrix(U, dimForce));
 482             fvVectorMatrix& fvm = tfvm();
 483
 484             fvm.source() = -UTrans()/(this->db().time().deltaT());
 485
 486             return tfvm;
 487         }
 488     }
 489
 490     return tmp<fvVectorMatrix>(new fvVectorMatrix(U, dimForce));
 491 }
 492
 493
 494 template<class CloudType>
 495 inline const Foam::tmp<Foam::volScalarField>
 496 Foam::KinematicCloud<CloudType>::theta() const
 497 {
 498     tmp<volScalarField> ttheta
 499     (
 500         new volScalarField
 501         (
 502             IOobject
 503             (
 504                 this->name() + "Theta",
 505                 this->db().time().timeName(),
 506                 this->db(),
 507                 IOobject::NO_READ,
 508                 IOobject::NO_WRITE,
 509                 false
 510             ),
 511             mesh_,
 512             dimensionedScalar("zero", dimless, 0.0),
 513             zeroGradientFvPatchScalarField::typeName
 514         )
 515     );
 516
 517     volScalarField& theta = ttheta();
 518     forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
 519     {
 520         const parcelType& p = iter();
 521         const label cellI = p.cell();
 522
 523         theta[cellI] += p.nParticle()*p.volume();
 524     }
 525
 526     theta.internalField() /= mesh_.V();
 527     theta.correctBoundaryConditions();
 528
 529     return ttheta;
 530 }
 531
 532
 533 template<class CloudType>
 534 inline const Foam::tmp<Foam::volScalarField>
 535 Foam::KinematicCloud<CloudType>::alpha() const
 536 {
 537     tmp<volScalarField> talpha
 538     (
 539         new volScalarField
 540         (
 541             IOobject
 542             (
 543                 this->name() + "Alpha",
 544                 this->db().time().timeName(),
 545                 this->db(),
 546                 IOobject::NO_READ,
 547                 IOobject::NO_WRITE,
 548                 false
 549             ),
 550             mesh_,
 551             dimensionedScalar("zero", dimless, 0.0)
 552         )
 553     );
 554
 555     scalarField& alpha = talpha().internalField();
 556     forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
 557     {
 558         const parcelType& p = iter();
 559         const label cellI = p.cell();
 560
 561         alpha[cellI] += p.nParticle()*p.mass();
 562     }
 563
 564     alpha /= (mesh_.V()*rho_);
 565
 566     return talpha;
 567 }
 568
 569
 570 template<class CloudType>
 571 inline const Foam::tmp<Foam::volScalarField>
 572 Foam::KinematicCloud<CloudType>::rhoEff() const
 573 {
 574     tmp<volScalarField> trhoEff
 575     (
 576         new volScalarField
 577         (
 578             IOobject
 579             (
 580                 this->name() + "RhoEff",
 581                 this->db().time().timeName(),
 582                 this->db(),
 583                 IOobject::NO_READ,
 584                 IOobject::NO_WRITE,
 585                 false
 586             ),
 587             mesh_,
 588             dimensionedScalar("zero", dimDensity, 0.0)
 589         )
 590     );
 591
 592     scalarField& rhoEff = trhoEff().internalField();
 593     forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
 594     {
 595         const parcelType& p = iter();
 596         const label cellI = p.cell();
 597
 598         rhoEff[cellI] += p.nParticle()*p.mass();
 599     }
 600
 601     rhoEff /= mesh_.V();
 602
 603     return trhoEff;
 604 }
 605
 606
 607 // ************************************************************************* //