multiphaseMixtureThermo.C 27.7 KB
Newer Older
1 2 3 4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
OpenFOAM bot's avatar
OpenFOAM bot committed
5
    \\  /    A nd           | www.openfoam.com
6
     \\/     M anipulation  |
OpenFOAM bot's avatar
OpenFOAM bot committed
7
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8
    Copyright (C) 2013-2017 OpenFOAM Foundation
9
    Copyright (C) 2020 OpenCFD Ltd.
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "multiphaseMixtureThermo.H"
#include "alphaContactAngleFvPatchScalarField.H"
#include "Time.H"
#include "subCycle.H"
#include "MULES.H"
#include "fvcDiv.H"
#include "fvcGrad.H"
#include "fvcSnGrad.H"
#include "fvcFlux.H"
#include "fvcMeshPhi.H"
#include "surfaceInterpolate.H"
40
#include "unitConversion.H"
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
    defineTypeNameAndDebug(multiphaseMixtureThermo, 0);
}


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

void Foam::multiphaseMixtureThermo::calcAlphas()
{
    scalar level = 0.0;
    alphas_ == 0.0;

57
    for (const phaseModel& phase : phases_)
58
    {
59
        alphas_ += level * phase;
60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83
        level += 1.0;
    }
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::multiphaseMixtureThermo::multiphaseMixtureThermo
(
    const volVectorField& U,
    const surfaceScalarField& phi
)
:
    psiThermo(U.mesh(), word::null),
    phases_(lookup("phases"), phaseModel::iNew(p_, T_)),

    mesh_(U.mesh()),
    U_(U),
    phi_(phi),

    rhoPhi_
    (
        IOobject
        (
84
            "rhoPhi",
85 86 87 88 89 90
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh_,
91
        dimensionedScalar(dimMass/dimTime, Zero)
92 93 94 95 96 97 98 99 100 101 102 103 104
    ),

    alphas_
    (
        IOobject
        (
            "alphas",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_,
105
        dimensionedScalar(dimless, Zero)
106 107 108 109 110 111 112
    ),

    sigmas_(lookup("sigmas")),
    dimSigma_(1, 0, -2, 0, 0),
    deltaN_
    (
        "deltaN",
113
        1e-8/cbrt(average(mesh_.V()))
114 115
    )
{
116
    rhoPhi_.setOriented();
117 118 119 120 121 122 123 124 125 126
    calcAlphas();
    alphas_.write();
    correct();
}


// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //

void Foam::multiphaseMixtureThermo::correct()
{
127
    for (phaseModel& phase : phases_)
128
    {
129
        phase.correct();
130 131
    }

132
    auto phasei = phases_.cbegin();
133 134 135 136 137

    psi_ = phasei()*phasei().thermo().psi();
    mu_ = phasei()*phasei().thermo().mu();
    alpha_ = phasei()*phasei().thermo().alpha();

138
    for (++phasei; phasei != phases_.cend(); ++phasei)
139 140 141 142 143 144 145 146 147 148
    {
        psi_ += phasei()*phasei().thermo().psi();
        mu_ += phasei()*phasei().thermo().mu();
        alpha_ += phasei()*phasei().thermo().alpha();
    }
}


void Foam::multiphaseMixtureThermo::correctRho(const volScalarField& dp)
{
149
    for (phaseModel& phase : phases_)
150
    {
151
        phase.thermo().rho() += phase.thermo().psi()*dp;
152 153 154 155
    }
}


156 157
Foam::word Foam::multiphaseMixtureThermo::thermoName() const
{
158
    auto phasei = phases_.cbegin();
159 160 161

    word name = phasei().thermo().thermoName();

162
    for (++phasei; phasei != phases_.cend(); ++phasei)
163 164 165 166 167 168 169 170
    {
        name += ',' + phasei().thermo().thermoName();
    }

    return name;
}


171 172
bool Foam::multiphaseMixtureThermo::incompressible() const
{
173
    for (const phaseModel& phase : phases_)
174
    {
175 176 177 178
        if (!phase.thermo().incompressible())
        {
            return false;
        }
179 180
    }

181
    return true;
182 183 184 185 186
}


bool Foam::multiphaseMixtureThermo::isochoric() const
{
187
    for (const phaseModel& phase : phases_)
188
    {
189 190 191 192
        if (!phase.thermo().isochoric())
        {
            return false;
        }
193 194
    }

195
    return true;
196 197 198 199 200 201 202 203 204
}


Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::he
(
    const volScalarField& p,
    const volScalarField& T
) const
{
205
    auto phasei = phases_.cbegin();
206 207 208

    tmp<volScalarField> the(phasei()*phasei().thermo().he(p, T));

209
    for (++phasei; phasei != phases_.cend(); ++phasei)
210
    {
211
        the.ref() += phasei()*phasei().thermo().he(p, T);
212 213 214 215 216 217 218 219 220 221 222 223 224
    }

    return the;
}


Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::he
(
    const scalarField& p,
    const scalarField& T,
    const labelList& cells
) const
{
225
    auto phasei = phases_.cbegin();
226 227 228 229 230 231

    tmp<scalarField> the
    (
        scalarField(phasei(), cells)*phasei().thermo().he(p, T, cells)
    );

232
    for (++phasei; phasei != phases_.cend(); ++phasei)
233
    {
234 235
        the.ref() +=
            scalarField(phasei(), cells)*phasei().thermo().he(p, T, cells);
236 237 238 239 240 241 242 243 244 245 246 247 248
    }

    return the;
}


Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::he
(
    const scalarField& p,
    const scalarField& T,
    const label patchi
) const
{
249
    auto phasei = phases_.cbegin();
250 251 252 253 254 255

    tmp<scalarField> the
    (
        phasei().boundaryField()[patchi]*phasei().thermo().he(p, T, patchi)
    );

256
    for (++phasei; phasei != phases_.cend(); ++phasei)
257
    {
258
        the.ref() +=
259 260 261 262 263 264 265 266 267
            phasei().boundaryField()[patchi]*phasei().thermo().he(p, T, patchi);
    }

    return the;
}


Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::hc() const
{
268
    auto phasei = phases_.cbegin();
269 270 271

    tmp<volScalarField> thc(phasei()*phasei().thermo().hc());

272
    for (++phasei; phasei != phases_.cend(); ++phasei)
273
    {
274
        thc.ref() += phasei()*phasei().thermo().hc();
275 276 277 278 279 280 281 282 283 284 285 286 287 288
    }

    return thc;
}


Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::THE
(
    const scalarField& h,
    const scalarField& p,
    const scalarField& T0,
    const labelList& cells
) const
{
289
    NotImplemented;
290 291 292 293 294 295 296 297 298 299 300 301
    return T0;
}


Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::THE
(
    const scalarField& h,
    const scalarField& p,
    const scalarField& T0,
    const label patchi
) const
{
302
    NotImplemented;
303 304 305 306 307 308
    return T0;
}


Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::rho() const
{
309
    auto phasei = phases_.cbegin();
310 311 312

    tmp<volScalarField> trho(phasei()*phasei().thermo().rho());

313
    for (++phasei; phasei != phases_.cend(); ++phasei)
314
    {
315
        trho.ref() += phasei()*phasei().thermo().rho();
316 317 318 319 320 321
    }

    return trho;
}


322 323 324 325 326
Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::rho
(
    const label patchi
) const
{
327
    auto phasei = phases_.cbegin();
328 329 330 331 332 333

    tmp<scalarField> trho
    (
        phasei().boundaryField()[patchi]*phasei().thermo().rho(patchi)
    );

334
    for (++phasei; phasei != phases_.cend(); ++phasei)
335
    {
336
        trho.ref() +=
337 338 339 340 341 342 343
            phasei().boundaryField()[patchi]*phasei().thermo().rho(patchi);
    }

    return trho;
}


344 345
Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::Cp() const
{
346
    auto phasei = phases_.cbegin();
347 348 349

    tmp<volScalarField> tCp(phasei()*phasei().thermo().Cp());

350
    for (++phasei; phasei != phases_.cend(); ++phasei)
351
    {
352
        tCp.ref() += phasei()*phasei().thermo().Cp();
353 354 355 356 357 358 359 360 361 362 363 364 365
    }

    return tCp;
}


Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::Cp
(
    const scalarField& p,
    const scalarField& T,
    const label patchi
) const
{
366
    auto phasei = phases_.cbegin();
367 368 369 370 371 372

    tmp<scalarField> tCp
    (
        phasei().boundaryField()[patchi]*phasei().thermo().Cp(p, T, patchi)
    );

373
    for (++phasei; phasei != phases_.cend(); ++phasei)
374
    {
375
        tCp.ref() +=
376 377 378 379 380 381 382 383 384
            phasei().boundaryField()[patchi]*phasei().thermo().Cp(p, T, patchi);
    }

    return tCp;
}


Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::Cv() const
{
385
    auto phasei = phases_.cbegin();
386 387 388

    tmp<volScalarField> tCv(phasei()*phasei().thermo().Cv());

389
    for (++phasei; phasei != phases_.cend(); ++phasei)
390
    {
391
        tCv.ref() += phasei()*phasei().thermo().Cv();
392 393 394 395 396 397 398 399 400 401 402 403 404
    }

    return tCv;
}


Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::Cv
(
    const scalarField& p,
    const scalarField& T,
    const label patchi
) const
{
405
    auto phasei = phases_.cbegin();
406 407 408 409 410 411

    tmp<scalarField> tCv
    (
        phasei().boundaryField()[patchi]*phasei().thermo().Cv(p, T, patchi)
    );

412
    for (++phasei; phasei != phases_.cend(); ++phasei)
413
    {
414
        tCv.ref() +=
415 416 417 418 419 420 421 422 423
            phasei().boundaryField()[patchi]*phasei().thermo().Cv(p, T, patchi);
    }

    return tCv;
}


Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::gamma() const
{
424
    auto phasei = phases_.cbegin();
425 426 427

    tmp<volScalarField> tgamma(phasei()*phasei().thermo().gamma());

428
    for (++phasei; phasei != phases_.cend(); ++phasei)
429
    {
430
        tgamma.ref() += phasei()*phasei().thermo().gamma();
431 432 433 434 435 436 437 438 439 440 441 442 443
    }

    return tgamma;
}


Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::gamma
(
    const scalarField& p,
    const scalarField& T,
    const label patchi
) const
{
444
    auto phasei = phases_.cbegin();
445 446 447 448 449 450

    tmp<scalarField> tgamma
    (
        phasei().boundaryField()[patchi]*phasei().thermo().gamma(p, T, patchi)
    );

451
    for (++phasei; phasei != phases_.cend(); ++phasei)
452
    {
453
        tgamma.ref() +=
454 455 456 457 458 459 460 461 462 463
            phasei().boundaryField()[patchi]
           *phasei().thermo().gamma(p, T, patchi);
    }

    return tgamma;
}


Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::Cpv() const
{
464
    auto phasei = phases_.cbegin();
465 466 467

    tmp<volScalarField> tCpv(phasei()*phasei().thermo().Cpv());

468
    for (++phasei; phasei != phases_.cend(); ++phasei)
469
    {
470
        tCpv.ref() += phasei()*phasei().thermo().Cpv();
471 472 473 474 475 476 477 478 479 480 481 482 483
    }

    return tCpv;
}


Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::Cpv
(
    const scalarField& p,
    const scalarField& T,
    const label patchi
) const
{
484
    auto phasei = phases_.cbegin();
485 486 487 488 489 490

    tmp<scalarField> tCpv
    (
        phasei().boundaryField()[patchi]*phasei().thermo().Cpv(p, T, patchi)
    );

491
    for (++phasei; phasei != phases_.cend(); ++phasei)
492
    {
493
        tCpv.ref() +=
494 495 496 497 498 499 500 501 502 503
            phasei().boundaryField()[patchi]
           *phasei().thermo().Cpv(p, T, patchi);
    }

    return tCpv;
}


Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::CpByCpv() const
{
504
    auto phasei = phases_.cbegin();
505 506 507

    tmp<volScalarField> tCpByCpv(phasei()*phasei().thermo().CpByCpv());

508
    for (++phasei; phasei != phases_.cend(); ++phasei)
509
    {
510
        tCpByCpv.ref() += phasei()*phasei().thermo().CpByCpv();
511 512 513 514 515 516 517 518 519 520 521 522 523
    }

    return tCpByCpv;
}


Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::CpByCpv
(
    const scalarField& p,
    const scalarField& T,
    const label patchi
) const
{
524
    auto phasei = phases_.cbegin();
525 526 527 528 529 530

    tmp<scalarField> tCpByCpv
    (
        phasei().boundaryField()[patchi]*phasei().thermo().CpByCpv(p, T, patchi)
    );

531
    for (++phasei; phasei != phases_.cend(); ++phasei)
532
    {
533
        tCpByCpv.ref() +=
534 535 536 537 538 539 540 541
            phasei().boundaryField()[patchi]
           *phasei().thermo().CpByCpv(p, T, patchi);
    }

    return tCpByCpv;
}


542 543
Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::W() const
{
544
    auto phasei = phases_.cbegin();
545 546 547

    tmp<volScalarField> tW(phasei()*phasei().thermo().W());

548
    for (++phasei; phasei != phases_.cend(); ++phasei)
549 550 551 552 553 554 555 556
    {
        tW.ref() += phasei()*phasei().thermo().W();
    }

    return tW;
}


557 558 559 560 561 562 563 564 565 566 567 568 569 570 571
Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::nu() const
{
    return mu()/rho();
}


Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::nu
(
    const label patchi
) const
{
    return mu(patchi)/rho(patchi);
}


572 573
Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::kappa() const
{
574
    auto phasei = phases_.cbegin();
575 576 577

    tmp<volScalarField> tkappa(phasei()*phasei().thermo().kappa());

578
    for (++phasei; phasei != phases_.cend(); ++phasei)
579
    {
580
        tkappa.ref() += phasei()*phasei().thermo().kappa();
581 582 583 584 585 586 587 588 589 590 591
    }

    return tkappa;
}


Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::kappa
(
    const label patchi
) const
{
592
    auto phasei = phases_.cbegin();
593 594 595 596 597 598

    tmp<scalarField> tkappa
    (
        phasei().boundaryField()[patchi]*phasei().thermo().kappa(patchi)
    );

599
    for (++phasei; phasei != phases_.cend(); ++phasei)
600
    {
601
        tkappa.ref() +=
602 603 604 605 606 607 608
            phasei().boundaryField()[patchi]*phasei().thermo().kappa(patchi);
    }

    return tkappa;
}


609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647
Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::alphahe() const
{
    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();

    tmp<volScalarField> talphaEff(phasei()*phasei().thermo().alphahe());

    for (++phasei; phasei != phases_.end(); ++phasei)
    {
        talphaEff.ref() += phasei()*phasei().thermo().alphahe();
    }

    return talphaEff;
}


Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::alphahe
(
    const label patchi
) const
{
    PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();

    tmp<scalarField> talphaEff
    (
        phasei().boundaryField()[patchi]
       *phasei().thermo().alphahe(patchi)
    );

    for (++phasei; phasei != phases_.end(); ++phasei)
    {
        talphaEff.ref() +=
            phasei().boundaryField()[patchi]
           *phasei().thermo().alphahe(patchi);
    }

    return talphaEff;
}


648 649 650 651 652
Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::kappaEff
(
    const volScalarField& alphat
) const
{
653
    auto phasei = phases_.cbegin();
654 655 656

    tmp<volScalarField> tkappaEff(phasei()*phasei().thermo().kappaEff(alphat));

657
    for (++phasei; phasei != phases_.cend(); ++phasei)
658
    {
659
        tkappaEff.ref() += phasei()*phasei().thermo().kappaEff(alphat);
660 661 662 663 664 665 666 667 668 669 670 671
    }

    return tkappaEff;
}


Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::kappaEff
(
    const scalarField& alphat,
    const label patchi
) const
{
672
    auto phasei = phases_.cbegin();
673 674 675 676 677 678 679

    tmp<scalarField> tkappaEff
    (
        phasei().boundaryField()[patchi]
       *phasei().thermo().kappaEff(alphat, patchi)
    );

680
    for (++phasei; phasei != phases_.cend(); ++phasei)
681
    {
682
        tkappaEff.ref() +=
683 684 685 686 687 688 689 690 691 692 693 694 695
            phasei().boundaryField()[patchi]
           *phasei().thermo().kappaEff(alphat, patchi);
    }

    return tkappaEff;
}


Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::alphaEff
(
    const volScalarField& alphat
) const
{
696
    auto phasei = phases_.cbegin();
697 698 699

    tmp<volScalarField> talphaEff(phasei()*phasei().thermo().alphaEff(alphat));

700
    for (++phasei; phasei != phases_.cend(); ++phasei)
701
    {
702
        talphaEff.ref() += phasei()*phasei().thermo().alphaEff(alphat);
703 704 705 706 707 708 709 710 711 712 713 714
    }

    return talphaEff;
}


Foam::tmp<Foam::scalarField> Foam::multiphaseMixtureThermo::alphaEff
(
    const scalarField& alphat,
    const label patchi
) const
{
715
    auto phasei = phases_.cbegin();
716 717 718 719 720 721 722

    tmp<scalarField> talphaEff
    (
        phasei().boundaryField()[patchi]
       *phasei().thermo().alphaEff(alphat, patchi)
    );

723
    for (++phasei; phasei != phases_.cend(); ++phasei)
724
    {
725
        talphaEff.ref() +=
726 727 728 729 730 731 732 733 734 735
            phasei().boundaryField()[patchi]
           *phasei().thermo().alphaEff(alphat, patchi);
    }

    return talphaEff;
}


Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::rCv() const
{
736
    auto phasei = phases_.cbegin();
737 738 739

    tmp<volScalarField> trCv(phasei()/phasei().thermo().Cv());

740
    for (++phasei; phasei != phases_.cend(); ++phasei)
741
    {
742
        trCv.ref() += phasei()/phasei().thermo().Cv();
743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762
    }

    return trCv;
}


Foam::tmp<Foam::surfaceScalarField>
Foam::multiphaseMixtureThermo::surfaceTensionForce() const
{
    tmp<surfaceScalarField> tstf
    (
        new surfaceScalarField
        (
            IOobject
            (
                "surfaceTensionForce",
                mesh_.time().timeName(),
                mesh_
            ),
            mesh_,
763
            dimensionedScalar(dimensionSet(1, -2, -2, 0, 0), Zero)
764 765 766
        )
    );

767
    surfaceScalarField& stf = tstf.ref();
768
    stf.setOriented();
769

770
    forAllConstIters(phases_, phase1)
771
    {
772
        const phaseModel& alpha1 = *phase1;
773

774
        auto phase2 = phase1;
775

776
        for (++phase2; phase2 != phases_.cend(); ++phase2)
777
        {
778
            const phaseModel& alpha2 = *phase2;
779

780
            auto sigma = sigmas_.cfind(interfacePair(alpha1, alpha2));
781

782
            if (!sigma.found())
783
            {
784 785
                FatalErrorInFunction
                    << "Cannot find interface " << interfacePair(alpha1, alpha2)
786 787 788 789
                    << " in list of sigma values"
                    << exit(FatalError);
            }

790
            stf += dimensionedScalar("sigma", dimSigma_, *sigma)
791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806
               *fvc::interpolate(K(alpha1, alpha2))*
                (
                    fvc::interpolate(alpha2)*fvc::snGrad(alpha1)
                  - fvc::interpolate(alpha1)*fvc::snGrad(alpha2)
                );
        }
    }

    return tstf;
}


void Foam::multiphaseMixtureThermo::solve()
{
    const Time& runTime = mesh_.time();

807
    const dictionary& alphaControls = mesh_.solverDict("alpha");
808 809
    label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
    scalar cAlpha(alphaControls.get<scalar>("cAlpha"));
810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883

    volScalarField& alpha = phases_.first();

    if (nAlphaSubCycles > 1)
    {
        surfaceScalarField rhoPhiSum(0.0*rhoPhi_);
        dimensionedScalar totalDeltaT = runTime.deltaT();

        for
        (
            subCycle<volScalarField> alphaSubCycle(alpha, nAlphaSubCycles);
            !(++alphaSubCycle).end();
        )
        {
            solveAlphas(cAlpha);
            rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi_;
        }

        rhoPhi_ = rhoPhiSum;
    }
    else
    {
        solveAlphas(cAlpha);
    }
}


Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixtureThermo::nHatfv
(
    const volScalarField& alpha1,
    const volScalarField& alpha2
) const
{
    /*
    // Cell gradient of alpha
    volVectorField gradAlpha =
        alpha2*fvc::grad(alpha1) - alpha1*fvc::grad(alpha2);

    // Interpolated face-gradient of alpha
    surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
    */

    surfaceVectorField gradAlphaf
    (
        fvc::interpolate(alpha2)*fvc::interpolate(fvc::grad(alpha1))
      - fvc::interpolate(alpha1)*fvc::interpolate(fvc::grad(alpha2))
    );

    // Face unit interface normal
    return gradAlphaf/(mag(gradAlphaf) + deltaN_);
}


Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixtureThermo::nHatf
(
    const volScalarField& alpha1,
    const volScalarField& alpha2
) const
{
    // Face unit interface normal flux
    return nHatfv(alpha1, alpha2) & mesh_.Sf();
}


// Correction for the boundary condition on the unit normal nHat on
// walls to produce the correct contact angle.

// The dynamic contact angle is calculated from the component of the
// velocity on the direction of the interface, parallel to the wall.

void Foam::multiphaseMixtureThermo::correctContactAngle
(
    const phaseModel& alpha1,
    const phaseModel& alpha2,
884
    surfaceVectorField::Boundary& nHatb
885 886
) const
{
887
    const volScalarField::Boundary& gbf
888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906
        = alpha1.boundaryField();

    const fvBoundaryMesh& boundary = mesh_.boundary();

    forAll(boundary, patchi)
    {
        if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
        {
            const alphaContactAngleFvPatchScalarField& acap =
                refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);

            vectorField& nHatPatch = nHatb[patchi];

            vectorField AfHatPatch
            (
                mesh_.Sf().boundaryField()[patchi]
               /mesh_.magSf().boundaryField()[patchi]
            );

907 908
            const auto tp =
                acap.thetaProps().cfind(interfacePair(alpha1, alpha2));
909

910
            if (!tp.found())
911
            {
912 913
                FatalErrorInFunction
                    << "Cannot find interface " << interfacePair(alpha1, alpha2)
914 915 916 917 918
                    << "\n    in table of theta properties for patch "
                    << acap.patch().name()
                    << exit(FatalError);
            }

919
            const bool matched = (tp.key().first() == alpha1.name());
920

921
            const scalar theta0 = degToRad(tp().theta0(matched));
922 923
            scalarField theta(boundary[patchi].size(), theta0);

924
            const scalar uTheta = tp().uTheta();
925 926 927 928

            // Calculate the dynamic contact angle if required
            if (uTheta > SMALL)
            {
929 930
                const scalar thetaA = degToRad(tp().thetaA(matched));
                const scalar thetaR = degToRad(tp().thetaR(matched));
931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990

                // Calculated the component of the velocity parallel to the wall
                vectorField Uwall
                (
                    U_.boundaryField()[patchi].patchInternalField()
                  - U_.boundaryField()[patchi]
                );
                Uwall -= (AfHatPatch & Uwall)*AfHatPatch;

                // Find the direction of the interface parallel to the wall
                vectorField nWall
                (
                    nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
                );

                // Normalise nWall
                nWall /= (mag(nWall) + SMALL);

                // Calculate Uwall resolved normal to the interface parallel to
                // the interface
                scalarField uwall(nWall & Uwall);

                theta += (thetaA - thetaR)*tanh(uwall/uTheta);
            }


            // Reset nHatPatch to correspond to the contact angle

            scalarField a12(nHatPatch & AfHatPatch);

            scalarField b1(cos(theta));

            scalarField b2(nHatPatch.size());

            forAll(b2, facei)
            {
                b2[facei] = cos(acos(a12[facei]) - theta[facei]);
            }

            scalarField det(1.0 - a12*a12);

            scalarField a((b1 - a12*b2)/det);
            scalarField b((b2 - a12*b1)/det);

            nHatPatch = a*AfHatPatch + b*nHatPatch;

            nHatPatch /= (mag(nHatPatch) + deltaN_.value());
        }
    }
}


Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::K
(
    const phaseModel& alpha1,
    const phaseModel& alpha2
) const
{
    tmp<surfaceVectorField> tnHatfv = nHatfv(alpha1, alpha2);

991
    correctContactAngle(alpha1, alpha2, tnHatfv.ref().boundaryFieldRef());
992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011

    // Simple expression for curvature
    return -fvc::div(tnHatfv & mesh_.Sf());
}


Foam::tmp<Foam::volScalarField>
Foam::multiphaseMixtureThermo::nearInterface() const
{
    tmp<volScalarField> tnearInt
    (
        new volScalarField
        (
            IOobject
            (
                "nearInterface",
                mesh_.time().timeName(),
                mesh_
            ),
            mesh_,
1012
            dimensionedScalar(dimless, Zero)
1013 1014 1015
        )
    );

1016
    for (const phaseModel& phase : phases_)
1017
    {
1018
        tnearInt.ref() =
1019
            max(tnearInt(), pos0(phase - 0.01)*pos0(0.99 - phase));
1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030
    }

    return tnearInt;
}


void Foam::multiphaseMixtureThermo::solveAlphas
(
    const scalar cAlpha
)
{
1031 1032
    static label nSolves(-1);
    ++nSolves;
1033

1034 1035
    const word alphaScheme("div(phi,alpha)");
    const word alpharScheme("div(phirb,alpha)");
1036 1037 1038 1039

    surfaceScalarField phic(mag(phi_/mesh_.magSf()));
    phic = min(cAlpha*phic, max(phic));

1040
    PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
1041

1042 1043
    int phasei = 0;
    for (phaseModel& alpha : phases_)
1044
    {
1045
        alphaPhiCorrs.set
1046 1047 1048 1049
        (
            phasei,
            new surfaceScalarField
            (
1050
                phi_.name() + alpha.name(),
1051 1052 1053 1054 1055 1056 1057 1058 1059
                fvc::flux
                (
                    phi_,
                    alpha,
                    alphaScheme
                )
            )
        );

1060
        surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
1061

1062
        for (phaseModel& alpha2 : phases_)
1063 1064 1065 1066 1067
        {
            if (&alpha2 == &alpha) continue;

            surfaceScalarField phir(phic*nHatf(alpha, alpha2));

1068
            alphaPhiCorr += fvc::flux
1069 1070 1071 1072 1073 1074 1075 1076 1077
            (
                -fvc::flux(-phir, alpha2, alpharScheme),
                alpha,
                alpharScheme
            );
        }

        MULES::limit
        (
Henry's avatar
Henry committed
1078
            1.0/mesh_.time().deltaT().value(),
1079 1080 1081
            geometricOneField(),
            alpha,
            phi_,
1082
            alphaPhiCorr,
1083 1084
            zeroField(),
            zeroField(),
1085 1086
            oneField(),
            zeroField(),
1087 1088 1089
            true
        );

1090
        ++phasei;
1091 1092
    }

1093
    MULES::limitSum(alphaPhiCorrs);
1094

1095
    rhoPhi_ = dimensionedScalar(dimensionSet(1, 0, -1, 0, 0), Zero);
1096 1097 1098 1099 1100 1101 1102 1103 1104 1105

    volScalarField sumAlpha
    (
        IOobject
        (
            "sumAlpha",
            mesh_.time().timeName(),
            mesh_
        ),
        mesh_,
1106
        dimensionedScalar(dimless, Zero)
1107 1108 1109 1110 1111 1112 1113 1114
    );


    volScalarField divU(fvc::div(fvc::absolute(phi_, U_)));


    phasei = 0;

1115
    for (phaseModel& alpha : phases_)
1116
    {
1117 1118
        surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
        alphaPhi += upwind<scalar>(mesh_, phi_).flux(alpha);
1119

1120
        volScalarField::Internal Sp
1121 1122 1123 1124 1125 1126 1127 1128
        (
            IOobject
            (
                "Sp",
                mesh_.time().timeName(),
                mesh_
            ),
            mesh_,