applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.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 "multiphaseMixture.H"
  27 #include "alphaContactAngleFvPatchScalarField.H"
  28 #include "Time.H"
  29 #include "subCycle.H"
  30 #include "fvCFD.H"
  31
  32 // * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //
  33
  34 const scalar Foam::multiphaseMixture::convertToRad =
  35     Foam::constant::mathematical::pi/180.0;
  36
  37
  38 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  39
  40 void Foam::multiphaseMixture::calcAlphas()
  41 {
  42     scalar level = 0.0;
  43     alphas_ == 0.0;
  44
  45     forAllIter(PtrDictionary<phase>, phases_, iter)
  46     {
  47         alphas_ += level*iter();
  48         level += 1.0;
  49     }
  50
  51     alphas_.correctBoundaryConditions();
  52 }
  53
  54
  55 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  56
  57 Foam::multiphaseMixture::multiphaseMixture
  58 (
  59     const volVectorField& U,
  60     const surfaceScalarField& phi
  61 )
  62 :
  63     transportModel(U, phi),
  64     phases_(lookup("phases"), phase::iNew(U, phi)),
  65     refPhase_(*phases_.lookup(word(lookup("refPhase")))),
  66
  67     mesh_(U.mesh()),
  68     U_(U),
  69     phi_(phi),
  70
  71     rhoPhi_
  72     (
  73         IOobject
  74         (
  75             "rho*phi",
  76             mesh_.time().timeName(),
  77             mesh_,
  78             IOobject::NO_READ,
  79             IOobject::NO_WRITE
  80         ),
  81         mesh_,
  82         dimensionedScalar("rho*phi", dimMass/dimTime, 0.0)
  83     ),
  84
  85     alphas_
  86     (
  87         IOobject
  88         (
  89             "alphas",
  90             mesh_.time().timeName(),
  91             mesh_,
  92             IOobject::NO_READ,
  93             IOobject::AUTO_WRITE
  94         ),
  95         mesh_,
  96         dimensionedScalar("alphas", dimless, 0.0),
  97         zeroGradientFvPatchScalarField::typeName
  98     ),
  99
 100     sigmas_(lookup("sigmas")),
 101     dimSigma_(1, 0, -2, 0, 0),
 102     deltaN_
 103     (
 104         "deltaN",
 105         1e-8/pow(average(mesh_.V()), 1.0/3.0)
 106     )
 107 {
 108     calcAlphas();
 109     alphas_.write();
 110
 111     forAllIter(PtrDictionary<phase>, phases_, iter)
 112     {
 113         alphaTable_.add(iter());
 114     }
 115 }
 116
 117
 118 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 119
 120 Foam::tmp<Foam::volScalarField> Foam::multiphaseMixture::rho() const
 121 {
 122     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 123
 124     tmp<volScalarField> trho = iter()*iter().rho();
 125
 126     for (++iter; iter != phases_.end(); ++iter)
 127     {
 128         trho() += iter()*iter().rho();
 129     }
 130
 131     return trho;
 132 }
 133
 134
 135 Foam::tmp<Foam::volScalarField> Foam::multiphaseMixture::mu() const
 136 {
 137     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 138
 139     tmp<volScalarField> tmu = iter()*iter().rho()*iter().nu();
 140
 141     for (++iter; iter != phases_.end(); ++iter)
 142     {
 143         tmu() += iter()*iter().rho()*iter().nu();
 144     }
 145
 146     return tmu;
 147 }
 148
 149
 150 Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixture::muf() const
 151 {
 152     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 153
 154     tmp<surfaceScalarField> tmuf =
 155         fvc::interpolate(iter())*iter().rho()*fvc::interpolate(iter().nu());
 156
 157     for (++iter; iter != phases_.end(); ++iter)
 158     {
 159         tmuf() +=
 160             fvc::interpolate(iter())*iter().rho()*fvc::interpolate(iter().nu());
 161     }
 162
 163     return tmuf;
 164 }
 165
 166
 167 Foam::tmp<Foam::volScalarField> Foam::multiphaseMixture::nu() const
 168 {
 169     return mu()/rho();
 170 }
 171
 172
 173 Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixture::nuf() const
 174 {
 175     return muf()/fvc::interpolate(rho());
 176 }
 177
 178
 179 Foam::tmp<Foam::surfaceScalarField>
 180 Foam::multiphaseMixture::surfaceTensionForce() const
 181 {
 182     tmp<surfaceScalarField> tstf
 183     (
 184         new surfaceScalarField
 185         (
 186             IOobject
 187             (
 188                 "surfaceTensionForce",
 189                 mesh_.time().timeName(),
 190                 mesh_
 191             ),
 192             mesh_,
 193             dimensionedScalar
 194             (
 195                 "surfaceTensionForce",
 196                 dimensionSet(1, -2, -2, 0, 0),
 197                 0.0
 198             )
 199         )
 200     );
 201
 202     surfaceScalarField& stf = tstf();
 203
 204     forAllConstIter(PtrDictionary<phase>, phases_, iter1)
 205     {
 206         const phase& alpha1 = iter1();
 207
 208         PtrDictionary<phase>::const_iterator iter2 = iter1;
 209         ++iter2;
 210
 211         for (; iter2 != phases_.end(); ++iter2)
 212         {
 213             const phase& alpha2 = iter2();
 214
 215             sigmaTable::const_iterator sigma =
 216                 sigmas_.find(interfacePair(alpha1, alpha2));
 217
 218             if (sigma == sigmas_.end())
 219             {
 220                 FatalErrorIn("multiphaseMixture::surfaceTensionForce() const")
 221                     << "Cannot find interface " << interfacePair(alpha1, alpha2)
 222                     << " in list of sigma values"
 223                     << exit(FatalError);
 224             }
 225
 226             stf += dimensionedScalar("sigma", dimSigma_, sigma())
 227                *fvc::interpolate(K(alpha1, alpha2))*
 228                 (
 229                     fvc::interpolate(alpha2)*fvc::snGrad(alpha1)
 230                   - fvc::interpolate(alpha1)*fvc::snGrad(alpha2)
 231                 );
 232         }
 233     }
 234
 235     return tstf;
 236 }
 237
 238
 239 void Foam::multiphaseMixture::solve()
 240 {
 241     forAllIter(PtrDictionary<phase>, phases_, iter)
 242     {
 243         iter().correct();
 244     }
 245
 246     const Time& runTime = mesh_.time();
 247
 248     const dictionary& pimpleDict = mesh_.solutionDict().subDict("PIMPLE");
 249
 250     label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
 251
 252     label nAlphaCorr(readLabel(pimpleDict.lookup("nAlphaCorr")));
 253
 254     bool cycleAlpha(Switch(pimpleDict.lookup("cycleAlpha")));
 255
 256     scalar cAlpha(readScalar(pimpleDict.lookup("cAlpha")));
 257
 258
 259     volScalarField& alpha = phases_.first();
 260
 261     if (nAlphaSubCycles > 1)
 262     {
 263         surfaceScalarField rhoPhiSum(0.0*rhoPhi_);
 264         dimensionedScalar totalDeltaT = runTime.deltaT();
 265
 266         for
 267         (
 268             subCycle<volScalarField> alphaSubCycle(alpha, nAlphaSubCycles);
 269             !(++alphaSubCycle).end();
 270         )
 271         {
 272             solveAlphas(nAlphaCorr, cycleAlpha, cAlpha);
 273             rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi_;
 274         }
 275
 276         rhoPhi_ = rhoPhiSum;
 277     }
 278     else
 279     {
 280         solveAlphas(nAlphaCorr, cycleAlpha, cAlpha);
 281     }
 282 }
 283
 284
 285 void Foam::multiphaseMixture::correct()
 286 {}
 287
 288
 289 Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixture::nHatfv
 290 (
 291     const volScalarField& alpha1,
 292     const volScalarField& alpha2
 293 ) const
 294 {
 295     /*
 296     // Cell gradient of alpha
 297     volVectorField gradAlpha =
 298         alpha2*fvc::grad(alpha1) - alpha1*fvc::grad(alpha2);
 299
 300     // Interpolated face-gradient of alpha
 301     surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
 302     */
 303
 304     surfaceVectorField gradAlphaf
 305     (
 306         fvc::interpolate(alpha2)*fvc::interpolate(fvc::grad(alpha1))
 307       - fvc::interpolate(alpha1)*fvc::interpolate(fvc::grad(alpha2))
 308     );
 309
 310     // Face unit interface normal
 311     return gradAlphaf/(mag(gradAlphaf) + deltaN_);
 312 }
 313
 314
 315 Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixture::nHatf
 316 (
 317     const volScalarField& alpha1,
 318     const volScalarField& alpha2
 319 ) const
 320 {
 321     // Face unit interface normal flux
 322     return nHatfv(alpha1, alpha2) & mesh_.Sf();
 323 }
 324
 325
 326 // Correction for the boundary condition on the unit normal nHat on
 327 // walls to produce the correct contact angle.
 328
 329 // The dynamic contact angle is calculated from the component of the
 330 // velocity on the direction of the interface, parallel to the wall.
 331
 332 void Foam::multiphaseMixture::correctContactAngle
 333 (
 334     const phase& alpha1,
 335     const phase& alpha2,
 336     surfaceVectorField::GeometricBoundaryField& nHatb
 337 ) const
 338 {
 339     const volScalarField::GeometricBoundaryField& gbf
 340         = refPhase_.boundaryField();
 341
 342     const fvBoundaryMesh& boundary = mesh_.boundary();
 343
 344     forAll(boundary, patchi)
 345     {
 346         if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
 347         {
 348             const alphaContactAngleFvPatchScalarField& acap =
 349                 refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);
 350
 351             vectorField& nHatPatch = nHatb[patchi];
 352
 353             vectorField AfHatPatch
 354             (
 355                 mesh_.Sf().boundaryField()[patchi]
 356                /mesh_.magSf().boundaryField()[patchi]
 357             );
 358
 359             alphaContactAngleFvPatchScalarField::thetaPropsTable::
 360                 const_iterator tp =
 361                 acap.thetaProps().find(interfacePair(alpha1, alpha2));
 362
 363             if (tp == acap.thetaProps().end())
 364             {
 365                 FatalErrorIn
 366                 (
 367                     "multiphaseMixture::correctContactAngle"
 368                     "(const phase& alpha1, const phase& alpha2, "
 369                     "fvPatchVectorFieldField& nHatb) const"
 370                 )   << "Cannot find interface " << interfacePair(alpha1, alpha2)
 371                     << "\n    in table of theta properties for patch "
 372                     << acap.patch().name()
 373                     << exit(FatalError);
 374             }
 375
 376             bool matched = (tp.key().first() == alpha1.name());
 377
 378             scalar theta0 = convertToRad*tp().theta0(matched);
 379             scalarField theta(boundary[patchi].size(), theta0);
 380
 381             scalar uTheta = tp().uTheta();
 382
 383             // Calculate the dynamic contact angle if required
 384             if (uTheta > SMALL)
 385             {
 386                 scalar thetaA = convertToRad*tp().thetaA(matched);
 387                 scalar thetaR = convertToRad*tp().thetaR(matched);
 388
 389                 // Calculated the component of the velocity parallel to the wall
 390                 vectorField Uwall
 391                 (
 392                     U_.boundaryField()[patchi].patchInternalField()
 393                   - U_.boundaryField()[patchi]
 394                 );
 395                 Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
 396
 397                 // Find the direction of the interface parallel to the wall
 398                 vectorField nWall
 399                 (
 400                     nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
 401                 );
 402
 403                 // Normalise nWall
 404                 nWall /= (mag(nWall) + SMALL);
 405
 406                 // Calculate Uwall resolved normal to the interface parallel to
 407                 // the interface
 408                 scalarField uwall(nWall & Uwall);
 409
 410                 theta += (thetaA - thetaR)*tanh(uwall/uTheta);
 411             }
 412
 413
 414             // Reset nHatPatch to correspond to the contact angle
 415
 416             scalarField a12(nHatPatch & AfHatPatch);
 417
 418             scalarField b1(cos(theta));
 419
 420             scalarField b2(nHatPatch.size());
 421
 422             forAll(b2, facei)
 423             {
 424                 b2[facei] = cos(acos(a12[facei]) - theta[facei]);
 425             }
 426
 427             scalarField det(1.0 - a12*a12);
 428
 429             scalarField a((b1 - a12*b2)/det);
 430             scalarField b((b2 - a12*b1)/det);
 431
 432             nHatPatch = a*AfHatPatch + b*nHatPatch;
 433
 434             nHatPatch /= (mag(nHatPatch) + deltaN_.value());
 435         }
 436     }
 437 }
 438
 439
 440 Foam::tmp<Foam::volScalarField> Foam::multiphaseMixture::K
 441 (
 442     const phase& alpha1,
 443     const phase& alpha2
 444 ) const
 445 {
 446     tmp<surfaceVectorField> tnHatfv = nHatfv(alpha1, alpha2);
 447
 448     correctContactAngle(alpha1, alpha2, tnHatfv().boundaryField());
 449
 450     // Simple expression for curvature
 451     return -fvc::div(tnHatfv & mesh_.Sf());
 452 }
 453
 454
 455 Foam::tmp<Foam::volScalarField>
 456 Foam::multiphaseMixture::nearInterface() const
 457 {
 458     tmp<volScalarField> tnearInt
 459     (
 460         new volScalarField
 461         (
 462             IOobject
 463             (
 464                 "nearInterface",
 465                 mesh_.time().timeName(),
 466                 mesh_
 467             ),
 468             mesh_,
 469             dimensionedScalar("nearInterface", dimless, 0.0)
 470         )
 471     );
 472
 473     forAllConstIter(PtrDictionary<phase>, phases_, iter)
 474     {
 475         tnearInt() = max(tnearInt(), pos(iter() - 0.01)*pos(0.99 - iter()));
 476     }
 477
 478     return tnearInt;
 479 }
 480
 481
 482 void Foam::multiphaseMixture::solveAlphas
 483 (
 484     const label nAlphaCorr,
 485     const bool cycleAlpha,
 486     const scalar cAlpha
 487 )
 488 {
 489     static label nSolves=-1;
 490     nSolves++;
 491
 492     word alphaScheme("div(phi,alpha)");
 493     word alphacScheme("div(phirb,alpha)");
 494
 495     tmp<fv::convectionScheme<scalar> > mvConvection
 496     (
 497         fv::convectionScheme<scalar>::New
 498         (
 499             mesh_,
 500             alphaTable_,
 501             phi_,
 502             mesh_.divScheme(alphaScheme)
 503         )
 504     );
 505
 506     surfaceScalarField phic(mag(phi_/mesh_.magSf()));
 507     phic = min(cAlpha*phic, max(phic));
 508
 509     for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
 510     {
 511         phase* refPhasePtr = &refPhase_;
 512
 513         if (cycleAlpha)
 514         {
 515             PtrDictionary<phase>::iterator refPhaseIter = phases_.begin();
 516             for (label i=0; i<nSolves%phases_.size(); i++)
 517             {
 518                 ++refPhaseIter;
 519             }
 520             refPhasePtr = &refPhaseIter();
 521         }
 522
 523         phase& refPhase = *refPhasePtr;
 524
 525         volScalarField refPhaseNew(refPhase);
 526         refPhaseNew == 1.0;
 527
 528         rhoPhi_ = phi_*refPhase.rho();
 529
 530         forAllIter(PtrDictionary<phase>, phases_, iter)
 531         {
 532             phase& alpha = iter();
 533
 534             if (&alpha == &refPhase) continue;
 535
 536             fvScalarMatrix alphaEqn
 537             (
 538                 fvm::ddt(alpha)
 539               + mvConvection->fvmDiv(phi_, alpha)
 540             );
 541
 542             forAllIter(PtrDictionary<phase>, phases_, iter2)
 543             {
 544                 phase& alpha2 = iter2();
 545
 546                 if (&alpha2 == &alpha) continue;
 547
 548                 surfaceScalarField phir(phic*nHatf(alpha, alpha2));
 549                 surfaceScalarField phirb12
 550                 (
 551                     -fvc::flux(-phir, alpha2, alphacScheme)
 552                 );
 553
 554                 alphaEqn += fvm::div(phirb12, alpha, alphacScheme);
 555             }
 556
 557             alphaEqn.solve(mesh_.solver("alpha"));
 558
 559             rhoPhi_ += alphaEqn.flux()*(alpha.rho() - refPhase.rho());
 560
 561             Info<< alpha.name() << " volume fraction, min, max = "
 562                 << alpha.weightedAverage(mesh_.V()).value()
 563                 << ' ' << min(alpha).value()
 564                 << ' ' << max(alpha).value()
 565                 << endl;
 566
 567             refPhaseNew == refPhaseNew - alpha;
 568         }
 569
 570         refPhase == refPhaseNew;
 571
 572         Info<< refPhase.name() << " volume fraction, min, max = "
 573             << refPhase.weightedAverage(mesh_.V()).value()
 574             << ' ' << min(refPhase).value()
 575             << ' ' << max(refPhase).value()
 576             << endl;
 577     }
 578
 579     calcAlphas();
 580 }
 581
 582
 583 bool Foam::multiphaseMixture::read()
 584 {
 585     if (transportModel::read())
 586     {
 587         bool readOK = true;
 588
 589         PtrList<entry> phaseData(lookup("phases"));
 590         label phasei = 0;
 591
 592         forAllIter(PtrDictionary<phase>, phases_, iter)
 593         {
 594             readOK &= iter().read(phaseData[phasei++].dict());
 595         }
 596
 597         lookup("sigmas") >> sigmas_;
 598
 599         return readOK;
 600     }
 601     else
 602     {
 603         return false;
 604     }
 605 }
 606
 607
 608 // ************************************************************************* //