src/regionModels/surfaceFilmModels/kinematicSingleLayer/kinematicSingleLayer.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 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 "kinematicSingleLayer.H"
  27 #include "fvm.H"
  28 #include "fvcDiv.H"
  29 #include "fvcLaplacian.H"
  30 #include "fvcSnGrad.H"
  31 #include "fvcReconstruct.H"
  32 #include "fvcVolumeIntegrate.H"
  33 #include "addToRunTimeSelectionTable.H"
  34 #include "directMappedWallPolyPatch.H"
  35 #include "mapDistribute.H"
  36
  37 #include "cachedRandom.H"
  38 #include "normal.H"
  39 #include "mathematicalConstants.H"
  40
  41 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  42
  43 namespace Foam
  44 {
  45 namespace regionModels
  46 {
  47 namespace surfaceFilmModels
  48 {
  49
  50 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  51
  52 defineTypeNameAndDebug(kinematicSingleLayer, 0);
  53
  54 addToRunTimeSelectionTable(surfaceFilmModel, kinematicSingleLayer, mesh);
  55
  56 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
  57
  58 bool kinematicSingleLayer::read()
  59 {
  60     if (surfaceFilmModel::read())
  61     {
  62         const dictionary& solution = this->solution().subDict("PISO");
  63         solution.lookup("momentumPredictor") >> momentumPredictor_;
  64         solution.lookup("nOuterCorr") >> nOuterCorr_;
  65         solution.lookup("nCorr") >> nCorr_;
  66         solution.lookup("nNonOrthCorr") >> nNonOrthCorr_;
  67
  68         return true;
  69     }
  70     else
  71     {
  72         return false;
  73     }
  74 }
  75
  76
  77 void kinematicSingleLayer::correctThermoFields()
  78 {
  79     if (thermoModel_ == tmConstant)
  80     {
  81         const dictionary& constDict(coeffs_.subDict("constantThermoCoeffs"));
  82         rho_ == dimensionedScalar(constDict.lookup("rho0"));
  83         mu_ == dimensionedScalar(constDict.lookup("mu0"));
  84         sigma_ == dimensionedScalar(constDict.lookup("sigma0"));
  85     }
  86     else
  87     {
  88         FatalErrorIn
  89         (
  90             "void Foam::surfaceFilmModels::kinematicSingleLayer::"
  91             "correctThermo()"
  92         )   << "Kinematic surface film must use "
  93             << thermoModelTypeNames_[thermoModel_] << "thermodynamics" << endl;
  94     }
  95 }
  96
  97
  98 void kinematicSingleLayer::resetPrimaryRegionSourceTerms()
  99 {
 100     if (debug)
 101     {
 102         Info<< "kinematicSingleLayer::resetPrimaryRegionSourceTerms()" << endl;
 103     }
 104
 105     rhoSpPrimary_ == dimensionedScalar("zero", rhoSp_.dimensions(), 0.0);
 106     USpPrimary_ == dimensionedVector("zero", USp_.dimensions(), vector::zero);
 107     pSpPrimary_ == dimensionedScalar("zero", pSp_.dimensions(), 0.0);
 108 }
 109
 110
 111 void kinematicSingleLayer::transferPrimaryRegionThermoFields()
 112 {
 113     if (debug)
 114     {
 115         Info<< "kinematicSingleLayer::"
 116             << "transferPrimaryRegionThermoFields()" << endl;
 117     }
 118     // Update fields from primary region via direct mapped
 119     // (coupled) boundary conditions
 120     UPrimary_.correctBoundaryConditions();
 121     pPrimary_.correctBoundaryConditions();
 122     rhoPrimary_.correctBoundaryConditions();
 123     muPrimary_.correctBoundaryConditions();
 124 }
 125
 126
 127 void kinematicSingleLayer::transferPrimaryRegionSourceFields()
 128 {
 129     if (debug)
 130     {
 131         Info<< "kinematicSingleLayer::"
 132             << "transferPrimaryRegionSourceFields()" << endl;
 133     }
 134
 135     // Retrieve the source fields from the primary region via direct mapped
 136     // (coupled) boundary conditions
 137     // - fields require transfer of values for both patch AND to push the
 138     //   values into the first layer of internal cells
 139     rhoSp_.correctBoundaryConditions();
 140     USp_.correctBoundaryConditions();
 141     pSp_.correctBoundaryConditions();
 142
 143     // Convert accummulated source terms into per unit area per unit time
 144     // Note: boundary values will still have original (neat) values
 145     const scalar deltaT = time_.deltaTValue();
 146     rhoSp_.field() /= magSf()*deltaT;
 147     USp_.field() /= magSf()*deltaT;
 148     pSp_.field() /= magSf()*deltaT;
 149 }
 150
 151
 152 tmp<volScalarField> kinematicSingleLayer::pu()
 153 {
 154     return tmp<volScalarField>
 155     (
 156         new volScalarField
 157         (
 158             IOobject
 159             (
 160                 "pu",
 161                 time_.timeName(),
 162                 regionMesh(),
 163                 IOobject::NO_READ,
 164                 IOobject::NO_WRITE
 165             ),
 166             pPrimary_                  // pressure (mapped from primary region)
 167           - pSp_                           // accumulated particle impingement
 168           - fvc::laplacian(sigma_, delta_) // surface tension
 169         )
 170     );
 171 }
 172
 173
 174 tmp<volScalarField> kinematicSingleLayer::pp()
 175 {
 176     return tmp<volScalarField>
 177     (
 178         new volScalarField
 179         (
 180             IOobject
 181             (
 182                 "pp",
 183                 time_.timeName(),
 184                 regionMesh(),
 185                 IOobject::NO_READ,
 186                 IOobject::NO_WRITE
 187             ),
 188            -rho_*gNormClipped() // hydrostatic effect only
 189         )
 190     );
 191 }
 192
 193
 194 void kinematicSingleLayer::updateSubmodels()
 195 {
 196     if (debug)
 197     {
 198         Info<< "kinematicSingleLayer::updateSubmodels()" << endl;
 199     }
 200
 201     // Update injection model - mass returned is mass available for injection
 202     injection_.correct(availableMass_, cloudMassTrans_, cloudDiameterTrans_);
 203
 204     // Update source fields
 205     const dimensionedScalar deltaT = time().deltaT();
 206     rhoSp_ += cloudMassTrans_/magSf()/deltaT;
 207 }
 208
 209
 210 void kinematicSingleLayer::continuityCheck()
 211 {
 212     const volScalarField deltaRho0(deltaRho_);
 213
 214     solveContinuity();
 215
 216     if (debug)
 217     {
 218         const volScalarField mass(deltaRho_*magSf());
 219         const dimensionedScalar totalMass =
 220             fvc::domainIntegrate(mass)
 221           + dimensionedScalar("SMALL", dimMass*dimVolume, ROOTVSMALL);
 222
 223         const scalar sumLocalContErr =
 224             (
 225                 fvc::domainIntegrate(mag(mass - magSf()*deltaRho0))/totalMass
 226             ).value();
 227
 228        const scalar globalContErr =
 229             (
 230                 fvc::domainIntegrate(mass - magSf()*deltaRho0)/totalMass
 231             ).value();
 232
 233         cumulativeContErr_ += globalContErr;
 234
 235         Info<< "Surface film: " << type() << nl
 236             << "    time step continuity errors: sum local = "
 237             << sumLocalContErr << ", global = " << globalContErr
 238             << ", cumulative = " << cumulativeContErr_ << endl;
 239     }
 240 }
 241
 242
 243 void kinematicSingleLayer::solveContinuity()
 244 {
 245     if (debug)
 246     {
 247         Info<< "kinematicSingleLayer::solveContinuity()" << endl;
 248     }
 249
 250     solve
 251     (
 252         fvm::ddt(deltaRho_)
 253       + fvc::div(phi_)
 254      ==
 255       - rhoSp_
 256     );
 257 }
 258
 259
 260 void kinematicSingleLayer::updateSurfaceVelocities()
 261 {
 262     // Push boundary film velocity values into internal field
 263     for (label i=0; i<intCoupledPatchIDs_.size(); i++)
 264     {
 265         label patchI = intCoupledPatchIDs_[i];
 266         const polyPatch& pp = regionMesh().boundaryMesh()[patchI];
 267         UIndirectList<vector>(Uw_, pp.faceCells()) =
 268             U_.boundaryField()[patchI];
 269     }
 270     Uw_ -= nHat()*(Uw_ & nHat());
 271     Uw_.correctBoundaryConditions();
 272
 273     // TODO: apply quadratic profile to determine surface velocity
 274     Us_ = U_;
 275     Us_.correctBoundaryConditions();
 276 }
 277
 278
 279 tmp<Foam::fvVectorMatrix> kinematicSingleLayer::solveMomentum
 280 (
 281     const volScalarField& pu,
 282     const volScalarField& pp
 283 )
 284 {
 285     if (debug)
 286     {
 287         Info<< "kinematicSingleLayer::solveMomentum()" << endl;
 288     }
 289
 290     updateSurfaceVelocities();
 291
 292     // Momentum
 293     tmp<fvVectorMatrix> tUEqn
 294     (
 295         fvm::ddt(deltaRho_, U_)
 296       + fvm::div(phi_, U_)
 297      ==
 298       - USp_
 299       - fvm::SuSp(rhoSp_, U_)
 300       + forces_.correct(U_)
 301     );
 302
 303     fvVectorMatrix& UEqn = tUEqn();
 304
 305     UEqn.relax();
 306
 307     if (momentumPredictor_)
 308     {
 309         solve
 310         (
 311             UEqn
 312          ==
 313             fvc::reconstruct
 314             (
 315               - fvc::interpolate(delta_)
 316               * (
 317                     regionMesh().magSf()
 318                   * (
 319                         fvc::snGrad(pu, "snGrad(p)")
 320                       + fvc::snGrad(pp, "snGrad(p)")*fvc::interpolate(delta_)
 321                       + fvc::snGrad(delta_)*fvc::interpolate(pp)
 322                     )
 323                   - (fvc::interpolate(rho_*gTan()) & regionMesh().Sf())
 324                 )
 325             )
 326         );
 327
 328         // Remove any patch-normal components of velocity
 329         U_ -= nHat()*(nHat() & U_);
 330         U_.correctBoundaryConditions();
 331     }
 332
 333     return tUEqn;
 334 }
 335
 336
 337 void kinematicSingleLayer::solveThickness
 338 (
 339     const volScalarField& pu,
 340     const volScalarField& pp,
 341     const fvVectorMatrix& UEqn
 342 )
 343 {
 344     if (debug)
 345     {
 346         Info<< "kinematicSingleLayer::solveThickness()" << endl;
 347     }
 348
 349     volScalarField rUA(1.0/UEqn.A());
 350     U_ = rUA*UEqn.H();
 351
 352     surfaceScalarField deltarUAf(fvc::interpolate(delta_*rUA));
 353     surfaceScalarField rhof(fvc::interpolate(rho_));
 354
 355     surfaceScalarField phiAdd
 356     (
 357         "phiAdd",
 358         regionMesh().magSf()
 359       * (
 360             fvc::snGrad(pu, "snGrad(p)")
 361           + fvc::snGrad(pp, "snGrad(p)")*fvc::interpolate(delta_)
 362         )
 363       - (fvc::interpolate(rho_*gTan()) & regionMesh().Sf())
 364     );
 365     constrainFilmField(phiAdd, 0.0);
 366
 367     surfaceScalarField phid
 368     (
 369         "phid",
 370         (fvc::interpolate(U_*rho_) & regionMesh().Sf())
 371       - deltarUAf*phiAdd*rhof
 372     );
 373     constrainFilmField(phid, 0.0);
 374
 375     surfaceScalarField ddrhorUAppf
 376     (
 377         "deltaCoeff",
 378         fvc::interpolate(delta_)*deltarUAf*rhof*fvc::interpolate(pp)
 379     );
 380 //    constrainFilmField(ddrhorUAppf, 0.0);
 381
 382     for (int nonOrth=0; nonOrth<=nNonOrthCorr_; nonOrth++)
 383     {
 384         // Film thickness equation
 385         fvScalarMatrix deltaEqn
 386         (
 387             fvm::ddt(rho_, delta_)
 388           + fvm::div(phid, delta_)
 389           - fvm::laplacian(ddrhorUAppf, delta_)
 390          ==
 391           - rhoSp_
 392         );
 393
 394         deltaEqn.solve();
 395
 396         if (nonOrth == nNonOrthCorr_)
 397         {
 398             phiAdd +=
 399                 fvc::interpolate(pp)
 400               * fvc::snGrad(delta_)
 401               * regionMesh().magSf();
 402
 403             phi_ == deltaEqn.flux();
 404         }
 405     }
 406
 407     // Bound film thickness by a minimum of zero
 408     delta_.max(0.0);
 409
 410     // Update U field
 411     U_ -= fvc::reconstruct(deltarUAf*phiAdd);
 412
 413     // Remove any patch-normal components of velocity
 414     U_ -= nHat()*(nHat() & U_);
 415
 416     U_.correctBoundaryConditions();
 417
 418     // Continuity check
 419     continuityCheck();
 420 }
 421
 422
 423 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 424
 425 kinematicSingleLayer::kinematicSingleLayer
 426 (
 427     const word& modelType,
 428     const fvMesh& mesh,
 429     const dimensionedVector& g,
 430     const bool readFields
 431 )
 432 :
 433     surfaceFilmModel(modelType, mesh, g),
 434
 435     momentumPredictor_(solution().subDict("PISO").lookup("momentumPredictor")),
 436     nOuterCorr_(readLabel(solution().subDict("PISO").lookup("nOuterCorr"))),
 437     nCorr_(readLabel(solution().subDict("PISO").lookup("nCorr"))),
 438     nNonOrthCorr_
 439     (
 440         readLabel(solution().subDict("PISO").lookup("nNonOrthCorr"))
 441     ),
 442
 443     cumulativeContErr_(0.0),
 444
 445     rho_
 446     (
 447         IOobject
 448         (
 449             "rhof",
 450             time().timeName(),
 451             regionMesh(),
 452             IOobject::NO_READ,
 453             IOobject::AUTO_WRITE
 454         ),
 455         regionMesh(),
 456         dimensionedScalar("zero", dimDensity, 0.0),
 457         zeroGradientFvPatchScalarField::typeName
 458     ),
 459     mu_
 460     (
 461         IOobject
 462         (
 463             "muf",
 464             time().timeName(),
 465             regionMesh(),
 466             IOobject::NO_READ,
 467             IOobject::AUTO_WRITE
 468         ),
 469         regionMesh(),
 470         dimensionedScalar("zero", dimPressure*dimTime, 0.0),
 471         zeroGradientFvPatchScalarField::typeName
 472     ),
 473     sigma_
 474     (
 475         IOobject
 476         (
 477             "sigmaf",
 478             time().timeName(),
 479             regionMesh(),
 480             IOobject::NO_READ,
 481             IOobject::AUTO_WRITE
 482         ),
 483         regionMesh(),
 484         dimensionedScalar("zero", dimMass/sqr(dimTime), 0.0),
 485         zeroGradientFvPatchScalarField::typeName
 486     ),
 487
 488     delta_
 489     (
 490         IOobject
 491         (
 492             "deltaf",
 493             time().timeName(),
 494             regionMesh(),
 495             IOobject::MUST_READ,
 496             IOobject::AUTO_WRITE
 497         ),
 498         regionMesh()
 499     ),
 500     U_
 501     (
 502         IOobject
 503         (
 504             "Uf",
 505             time().timeName(),
 506             regionMesh(),
 507             IOobject::MUST_READ,
 508             IOobject::AUTO_WRITE
 509         ),
 510         regionMesh()
 511     ),
 512     Us_
 513     (
 514         IOobject
 515         (
 516             "Usf",
 517             time().timeName(),
 518             regionMesh(),
 519             IOobject::NO_READ,
 520             IOobject::NO_WRITE
 521         ),
 522         U_,
 523         zeroGradientFvPatchScalarField::typeName
 524     ),
 525     Uw_
 526     (
 527         IOobject
 528         (
 529             "Uwf",
 530             time().timeName(),
 531             regionMesh(),
 532             IOobject::NO_READ,
 533             IOobject::NO_WRITE
 534         ),
 535         U_,
 536         zeroGradientFvPatchScalarField::typeName
 537     ),
 538     deltaRho_
 539     (
 540         IOobject
 541         (
 542             delta_.name() + "*" + rho_.name(),
 543             time().timeName(),
 544             regionMesh(),
 545             IOobject::NO_READ,
 546             IOobject::NO_WRITE
 547         ),
 548         regionMesh(),
 549         dimensionedScalar("zero", delta_.dimensions()*rho_.dimensions(), 0.0),
 550         zeroGradientFvPatchScalarField::typeName
 551     ),
 552
 553     phi_
 554     (
 555         IOobject
 556         (
 557             "phi",
 558             time().timeName(),
 559             regionMesh(),
 560             IOobject::READ_IF_PRESENT,
 561             IOobject::AUTO_WRITE
 562         ),
 563         regionMesh(),
 564         dimLength*dimMass/dimTime
 565     ),
 566
 567     primaryMassTrans_
 568     (
 569         IOobject
 570         (
 571             "primaryMassTrans",
 572             time().timeName(),
 573             regionMesh(),
 574             IOobject::NO_READ,
 575             IOobject::NO_WRITE
 576         ),
 577         regionMesh(),
 578         dimensionedScalar("zero", dimMass, 0.0),
 579         zeroGradientFvPatchScalarField::typeName
 580     ),
 581     cloudMassTrans_
 582     (
 583         IOobject
 584         (
 585             "cloudMassTrans",
 586             time().timeName(),
 587             regionMesh(),
 588             IOobject::NO_READ,
 589             IOobject::NO_WRITE
 590         ),
 591         regionMesh(),
 592         dimensionedScalar("zero", dimMass, 0.0),
 593         zeroGradientFvPatchScalarField::typeName
 594     ),
 595     cloudDiameterTrans_
 596     (
 597         IOobject
 598         (
 599             "cloudDiameterTrans",
 600             time().timeName(),
 601             regionMesh(),
 602             IOobject::NO_READ,
 603             IOobject::NO_WRITE
 604         ),
 605         regionMesh(),
 606         dimensionedScalar("zero", dimLength, -1.0),
 607         zeroGradientFvPatchScalarField::typeName
 608     ),
 609
 610     USp_
 611     (
 612         IOobject
 613         (
 614             "USpf",
 615             time().timeName(),
 616             regionMesh(),
 617             IOobject::NO_READ,
 618             IOobject::NO_WRITE
 619         ),
 620         regionMesh(),
 621         dimensionedVector
 622         (
 623             "zero", dimMass*dimVelocity/dimArea/dimTime, vector::zero
 624         ),
 625         this->mappedPushedFieldPatchTypes<vector>()
 626     ),
 627     pSp_
 628     (
 629         IOobject
 630         (
 631             "pSpf",
 632             time_.timeName(),
 633             regionMesh(),
 634             IOobject::NO_READ,
 635             IOobject::NO_WRITE
 636         ),
 637         regionMesh(),
 638         dimensionedScalar("zero", dimPressure, 0.0),
 639         this->mappedPushedFieldPatchTypes<scalar>()
 640     ),
 641     rhoSp_
 642     (
 643         IOobject
 644         (
 645             "rhoSpf",
 646             time_.timeName(),
 647             regionMesh(),
 648             IOobject::NO_READ,
 649             IOobject::NO_WRITE
 650         ),
 651         regionMesh(),
 652         dimensionedScalar("zero", dimMass/dimTime/dimArea, 0.0),
 653         this->mappedPushedFieldPatchTypes<scalar>()
 654     ),
 655
 656     USpPrimary_
 657     (
 658         IOobject
 659         (
 660             USp_.name(), // must have same name as USp_ to enable mapping
 661             time().timeName(),
 662             primaryMesh(),
 663             IOobject::NO_READ,
 664             IOobject::NO_WRITE
 665         ),
 666         primaryMesh(),
 667         dimensionedVector("zero", USp_.dimensions(), vector::zero)
 668     ),
 669     pSpPrimary_
 670     (
 671         IOobject
 672         (
 673             pSp_.name(), // must have same name as pSp_ to enable mapping
 674             time().timeName(),
 675             primaryMesh(),
 676             IOobject::NO_READ,
 677             IOobject::NO_WRITE
 678         ),
 679         primaryMesh(),
 680         dimensionedScalar("zero", pSp_.dimensions(), 0.0)
 681     ),
 682     rhoSpPrimary_
 683     (
 684         IOobject
 685         (
 686             rhoSp_.name(), // must have same name as rhoSp_ to enable mapping
 687             time().timeName(),
 688             primaryMesh(),
 689             IOobject::NO_READ,
 690             IOobject::NO_WRITE
 691         ),
 692         primaryMesh(),
 693         dimensionedScalar("zero", rhoSp_.dimensions(), 0.0)
 694     ),
 695
 696     UPrimary_
 697     (
 698         IOobject
 699         (
 700             "U", // must have same name as U to enable mapping
 701             time().timeName(),
 702             regionMesh(),
 703             IOobject::NO_READ,
 704             IOobject::NO_WRITE
 705         ),
 706         regionMesh(),
 707         dimensionedVector("zero", dimVelocity, vector::zero),
 708         this->mappedFieldAndInternalPatchTypes<vector>()
 709     ),
 710     pPrimary_
 711     (
 712         IOobject
 713         (
 714             "p", // must have same name as p to enable mapping
 715             time().timeName(),
 716             regionMesh(),
 717             IOobject::NO_READ,
 718             IOobject::NO_WRITE
 719         ),
 720         regionMesh(),
 721         dimensionedScalar("zero", dimPressure, 0.0),
 722         this->mappedFieldAndInternalPatchTypes<scalar>()
 723     ),
 724     rhoPrimary_
 725     (
 726         IOobject
 727         (
 728             "rho", // must have same name as rho to enable mapping
 729             time().timeName(),
 730             regionMesh(),
 731             IOobject::NO_READ,
 732             IOobject::NO_WRITE
 733         ),
 734         regionMesh(),
 735         dimensionedScalar("zero", dimDensity, 0.0),
 736         this->mappedFieldAndInternalPatchTypes<scalar>()
 737     ),
 738     muPrimary_
 739     (
 740         IOobject
 741         (
 742             "mu", // must have same name as mu to enable mapping
 743             time().timeName(),
 744             regionMesh(),
 745             IOobject::NO_READ,
 746             IOobject::NO_WRITE
 747         ),
 748         regionMesh(),
 749         dimensionedScalar("zero", dimPressure*dimTime, 0.0),
 750         this->mappedFieldAndInternalPatchTypes<scalar>()
 751     ),
 752
 753     availableMass_(regionMesh().nCells(), 0.0),
 754
 755     injection_(*this, coeffs_),
 756
 757     forces_(*this, coeffs_),
 758
 759     addedMassTotal_(0.0)
 760 {
 761     if (readFields)
 762     {
 763         transferPrimaryRegionThermoFields();
 764
 765         correctThermoFields();
 766
 767         deltaRho_ == delta_*rho_;
 768         phi_ = fvc::interpolate(deltaRho_*U_) & regionMesh().Sf();
 769     }
 770 }
 771
 772
 773 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 774
 775 kinematicSingleLayer::~kinematicSingleLayer()
 776 {}
 777
 778
 779 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 780
 781 void kinematicSingleLayer::addSources
 782 (
 783     const label patchI,
 784     const label faceI,
 785     const scalar massSource,
 786     const vector& momentumSource,
 787     const scalar pressureSource,
 788     const scalar energySource
 789 )
 790 {
 791     if (debug)
 792     {
 793         Info<< "\nSurface film: " << type() << ": adding to film source:" << nl
 794             << "    mass     = " << massSource << nl
 795             << "    momentum = " << momentumSource << nl
 796             << "    pressure = " << pressureSource << endl;
 797     }
 798
 799     rhoSpPrimary_.boundaryField()[patchI][faceI] -= massSource;
 800     USpPrimary_.boundaryField()[patchI][faceI] -= momentumSource;
 801     pSpPrimary_.boundaryField()[patchI][faceI] -= pressureSource;
 802
 803     addedMassTotal_ += massSource;
 804 }
 805
 806
 807 void kinematicSingleLayer::preEvolveRegion()
 808 {
 809     if (debug)
 810     {
 811         Info<< "kinematicSingleLayer::preEvolveRegion()" << endl;
 812     }
 813
 814     transferPrimaryRegionThermoFields();
 815
 816     correctThermoFields();
 817
 818     transferPrimaryRegionSourceFields();
 819
 820     // Reset transfer fields
 821 //    availableMass_ = mass();
 822     availableMass_ = netMass();
 823     cloudMassTrans_ == dimensionedScalar("zero", dimMass, 0.0);
 824     cloudDiameterTrans_ == dimensionedScalar("zero", dimLength, -1.0);
 825 }
 826
 827
 828 void kinematicSingleLayer::evolveRegion()
 829 {
 830     if (debug)
 831     {
 832         Info<< "kinematicSingleLayer::evolveRegion()" << endl;
 833     }
 834
 835     updateSubmodels();
 836
 837     // Solve continuity for deltaRho_
 838     solveContinuity();
 839
 840     // Implicit pressure source coefficient - constant
 841     tmp<volScalarField> tpp(this->pp());
 842
 843     for (int oCorr=0; oCorr<nOuterCorr_; oCorr++)
 844     {
 845         // Explicit pressure source contribution - varies with delta_
 846         tmp<volScalarField> tpu(this->pu());
 847
 848         // Solve for momentum for U_
 849         tmp<fvVectorMatrix> UEqn = solveMomentum(tpu(), tpp());
 850
 851         // Film thickness correction loop
 852         for (int corr=1; corr<=nCorr_; corr++)
 853         {
 854             // Solve thickness for delta_
 855             solveThickness(tpu(), tpp(), UEqn());
 856         }
 857     }
 858
 859     // Update deltaRho_ with new delta_
 860     deltaRho_ == delta_*rho_;
 861
 862     // Update film wall and surface velocities
 863     updateSurfaceVelocities();
 864
 865     // Reset source terms for next time integration
 866     resetPrimaryRegionSourceTerms();
 867 }
 868
 869
 870 scalar kinematicSingleLayer::CourantNumber() const
 871 {
 872     scalar CoNum = 0.0;
 873
 874     if (regionMesh().nInternalFaces() > 0)
 875     {
 876         const scalar deltaT = time_.deltaTValue();
 877
 878         const surfaceScalarField SfUfbyDelta
 879         (
 880             regionMesh().surfaceInterpolation::deltaCoeffs()*mag(phi_)
 881         );
 882         const surfaceScalarField rhoDelta(fvc::interpolate(rho_*delta_));
 883         const surfaceScalarField& magSf = regionMesh().magSf();
 884
 885         forAll(rhoDelta, i)
 886         {
 887             if (rhoDelta[i] > ROOTVSMALL)
 888             {
 889                 CoNum = max(CoNum, SfUfbyDelta[i]/rhoDelta[i]/magSf[i]*deltaT);
 890             }
 891         }
 892     }
 893
 894     reduce(CoNum, maxOp<scalar>());
 895
 896     Info<< "Film max Courant number: " << CoNum << endl;
 897
 898     return CoNum;
 899 }
 900
 901
 902 const volVectorField& kinematicSingleLayer::U() const
 903 {
 904     return U_;
 905 }
 906
 907
 908 const volVectorField& kinematicSingleLayer::Us() const
 909 {
 910     return Us_;
 911 }
 912
 913
 914 const volVectorField& kinematicSingleLayer::Uw() const
 915 {
 916     return Uw_;
 917 }
 918
 919
 920 const surfaceScalarField& kinematicSingleLayer::phi() const
 921 {
 922     return phi_;
 923 }
 924
 925
 926 const volScalarField& kinematicSingleLayer::rho() const
 927 {
 928     return rho_;
 929 }
 930
 931
 932 const volScalarField& kinematicSingleLayer::T() const
 933 {
 934     FatalErrorIn
 935     (
 936         "const volScalarField& kinematicSingleLayer::T() const"
 937     )   << "T field not available for " << type() << abort(FatalError);
 938
 939     return volScalarField::null();
 940 }
 941
 942
 943 const volScalarField& kinematicSingleLayer::Ts() const
 944 {
 945     FatalErrorIn
 946     (
 947         "const volScalarField& kinematicSingleLayer::Ts() const"
 948     )   << "Ts field not available for " << type() << abort(FatalError);
 949
 950     return volScalarField::null();
 951 }
 952
 953
 954 const volScalarField& kinematicSingleLayer::Tw() const
 955 {
 956     FatalErrorIn
 957     (
 958         "const volScalarField& kinematicSingleLayer::Tw() const"
 959     )   << "Tw field not available for " << type() << abort(FatalError);
 960
 961     return volScalarField::null();
 962 }
 963
 964
 965 const volScalarField& kinematicSingleLayer::Cp() const
 966 {
 967     FatalErrorIn
 968     (
 969         "const volScalarField& kinematicSingleLayer::Cp() const"
 970     )   << "Cp field not available for " << type() << abort(FatalError);
 971
 972     return volScalarField::null();
 973 }
 974
 975
 976 const volScalarField& kinematicSingleLayer::kappa() const
 977 {
 978     FatalErrorIn
 979     (
 980         "const volScalarField& kinematicSingleLayer::kappa() const"
 981     )   << "kappa field not available for " << type() << abort(FatalError);
 982
 983     return volScalarField::null();
 984 }
 985
 986
 987 tmp<volScalarField> kinematicSingleLayer::primaryMassTrans() const
 988 {
 989     return tmp<volScalarField>
 990     (
 991         new volScalarField
 992         (
 993             IOobject
 994             (
 995                 "kinematicSingleLayer::primaryMassTrans",
 996                 time().timeName(),
 997                 primaryMesh(),
 998                 IOobject::NO_READ,
 999                 IOobject::NO_WRITE,
1000                 false
1001             ),
1002             primaryMesh(),
1003             dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
1004         )
1005     );
1006 }
1007
1008
1009 const volScalarField& kinematicSingleLayer::cloudMassTrans() const
1010 {
1011     return cloudMassTrans_;
1012 }
1013
1014
1015 const volScalarField& kinematicSingleLayer::cloudDiameterTrans() const
1016 {
1017     return cloudDiameterTrans_;
1018 }
1019
1020
1021 void kinematicSingleLayer::info() const
1022 {
1023     Info<< "\nSurface film: " << type() << endl;
1024
1025     Info<< indent << "added mass         = "
1026         << returnReduce<scalar>(addedMassTotal_, sumOp<scalar>()) << nl
1027         << indent << "current mass       = "
1028         << gSum((deltaRho_*magSf())()) << nl
1029         << indent << "min/max(mag(U))    = " << min(mag(U_)).value() << ", "
1030         << max(mag(U_)).value() << nl
1031         << indent << "min/max(delta)     = " << min(delta_).value() << ", "
1032         << max(delta_).value() << nl;
1033
1034     injection_.info(Info);
1035 }
1036
1037
1038 tmp<DimensionedField<scalar, volMesh> > kinematicSingleLayer::Srho() const
1039 {
1040     return tmp<DimensionedField<scalar, volMesh> >
1041     (
1042         new DimensionedField<scalar, volMesh>
1043         (
1044             IOobject
1045             (
1046                 "kinematicSingleLayer::Srho",
1047                 time().timeName(),
1048                 primaryMesh(),
1049                 IOobject::NO_READ,
1050                 IOobject::NO_WRITE,
1051                 false
1052             ),
1053             primaryMesh(),
1054             dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
1055         )
1056     );
1057 }
1058
1059
1060 tmp<DimensionedField<scalar, volMesh> > kinematicSingleLayer::Srho
1061 (
1062     const label i
1063 ) const
1064 {
1065     return tmp<DimensionedField<scalar, volMesh> >
1066     (
1067         new DimensionedField<scalar, volMesh>
1068         (
1069             IOobject
1070             (
1071                 "kinematicSingleLayer::Srho(" + Foam::name(i) + ")",
1072                 time().timeName(),
1073                 primaryMesh(),
1074                 IOobject::NO_READ,
1075                 IOobject::NO_WRITE,
1076                 false
1077             ),
1078             primaryMesh(),
1079             dimensionedScalar("zero", dimMass/dimVolume/dimTime, 0.0)
1080         )
1081     );
1082 }
1083
1084
1085 tmp<DimensionedField<scalar, volMesh> > kinematicSingleLayer::Sh() const
1086 {
1087     return tmp<DimensionedField<scalar, volMesh> >
1088     (
1089         new DimensionedField<scalar, volMesh>
1090         (
1091             IOobject
1092             (
1093                 "kinematicSingleLayer::Sh",
1094                 time().timeName(),
1095                 primaryMesh(),
1096                 IOobject::NO_READ,
1097                 IOobject::NO_WRITE,
1098                 false
1099             ),
1100             primaryMesh(),
1101             dimensionedScalar("zero", dimEnergy/dimVolume/dimTime, 0.0)
1102         )
1103     );
1104 }
1105
1106
1107 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1108
1109 } // End namespace surfaceFilmModels
1110 } // End namespace regionModels
1111 } // End namespace Foam
1112
1113 // ************************************************************************* //