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