applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C

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