multiphaseSystem.C 22.4 KB
Newer Older
Henry's avatar
Henry committed
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
Henry's avatar
Henry committed
6
     \\/     M anipulation  |
OpenFOAM bot's avatar
OpenFOAM bot committed
7
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8
    Copyright (C) 2011-2018 OpenFOAM Foundation
9
    Copyright (C) 2020 OpenCFD Ltd.
Henry's avatar
Henry committed
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
-------------------------------------------------------------------------------
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 "multiphaseSystem.H"
#include "alphaContactAngleFvPatchScalarField.H"
31
#include "fixedValueFvsPatchFields.H"
Henry's avatar
Henry committed
32 33 34
#include "Time.H"
#include "subCycle.H"
#include "MULES.H"
35 36
#include "surfaceInterpolate.H"
#include "fvcGrad.H"
Henry's avatar
Henry committed
37
#include "fvcSnGrad.H"
38
#include "fvcDiv.H"
Henry's avatar
Henry committed
39
#include "fvcFlux.H"
40
#include "fvcAverage.H"
41

42
#include "unitConversion.H"
Henry's avatar
Henry committed
43 44 45 46 47 48 49 50

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

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

51
    for (const phaseModel& phase : phases_)
Henry's avatar
Henry committed
52
    {
53
        alphas_ += level * phase;
Henry's avatar
Henry committed
54 55 56 57 58 59 60
        level += 1.0;
    }
}


void Foam::multiphaseSystem::solveAlphas()
{
61
    PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
Henry's avatar
Henry committed
62 63
    int phasei = 0;

64
    for (phaseModel& phase : phases_)
Henry's avatar
Henry committed
65
    {
66
        volScalarField& alpha1 = phase;
Henry's avatar
Henry committed
67

68
        alphaPhiCorrs.set
Henry's avatar
Henry committed
69 70 71 72
        (
            phasei,
            new surfaceScalarField
            (
73
                "phi" + alpha1.name() + "Corr",
Henry's avatar
Henry committed
74 75 76
                fvc::flux
                (
                    phi_,
77
                    phase,
andy's avatar
andy committed
78
                    "div(phi," + alpha1.name() + ')'
Henry's avatar
Henry committed
79 80 81 82
                )
            )
        );

83
        surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
Henry's avatar
Henry committed
84

85
        for (phaseModel& phase2 : phases_)
Henry's avatar
Henry committed
86
        {
andy's avatar
andy committed
87
            volScalarField& alpha2 = phase2;
Henry's avatar
Henry committed
88

89
            if (&phase2 == &phase) continue;
Henry's avatar
Henry committed
90

91
            surfaceScalarField phir(phase.phi() - phase2.phi());
92

93
            const auto cAlpha = cAlphas_.cfind(interfacePair(phase, phase2));
Henry's avatar
Henry committed
94

95
            if (cAlpha.found())
96 97 98
            {
                surfaceScalarField phic
                (
99
                    (mag(phi_) + mag(phir))/mesh_.magSf()
100 101
                );

102
                phir += min(cAlpha()*phic, max(phic))*nHatf(phase, phase2);
103 104
            }

105
            const word phirScheme
andy's avatar
andy committed
106 107 108 109
            (
                "div(phir," + alpha2.name() + ',' + alpha1.name() + ')'
            );

110
            alphaPhiCorr += fvc::flux
Henry's avatar
Henry committed
111
            (
andy's avatar
andy committed
112
                -fvc::flux(-phir, phase2, phirScheme),
113
                phase,
andy's avatar
andy committed
114
                phirScheme
Henry's avatar
Henry committed
115 116 117
            );
        }

118
        phase.correctInflowOutflow(alphaPhiCorr);
119

Henry's avatar
Henry committed
120 121
        MULES::limit
        (
Henry's avatar
Henry committed
122
            1.0/mesh_.time().deltaT().value(),
Henry's avatar
Henry committed
123
            geometricOneField(),
124
            phase,
Henry's avatar
Henry committed
125
            phi_,
126
            alphaPhiCorr,
Henry's avatar
Henry committed
127 128
            zeroField(),
            zeroField(),
129 130
            oneField(),
            zeroField(),
Henry's avatar
Henry committed
131 132 133
            true
        );

134
        ++phasei;
Henry's avatar
Henry committed
135 136
    }

137
    MULES::limitSum(alphaPhiCorrs);
Henry's avatar
Henry committed
138 139 140 141 142 143 144 145 146 147

    volScalarField sumAlpha
    (
        IOobject
        (
            "sumAlpha",
            mesh_.time().timeName(),
            mesh_
        ),
        mesh_,
148
        dimensionedScalar("sumAlpha", dimless, 0)
Henry's avatar
Henry committed
149 150 151 152
    );

    phasei = 0;

153
    for (phaseModel& phase : phases_)
Henry's avatar
Henry committed
154
    {
155
        surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
156 157
        alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase);
        phase.correctInflowOutflow(alphaPhi);
Henry's avatar
Henry committed
158 159 160 161

        MULES::explicitSolve
        (
            geometricOneField(),
162
            phase,
163
            alphaPhi
Henry's avatar
Henry committed
164 165
        );

166
        phase.alphaPhi() = alphaPhi;
Henry's avatar
Henry committed
167

168 169 170 171
        Info<< phase.name() << " volume fraction, min, max = "
            << phase.weightedAverage(mesh_.V()).value()
            << ' ' << min(phase).value()
            << ' ' << max(phase).value()
Henry's avatar
Henry committed
172 173
            << endl;

174
        sumAlpha += phase;
Henry's avatar
Henry committed
175

176
        ++phasei;
Henry's avatar
Henry committed
177 178 179 180 181 182 183 184
    }

    Info<< "Phase-sum volume fraction, min, max = "
        << sumAlpha.weightedAverage(mesh_.V()).value()
        << ' ' << min(sumAlpha).value()
        << ' ' << max(sumAlpha).value()
        << endl;

185 186
    // Correct the sum of the phase-fractions to avoid 'drift'
    volScalarField sumCorr(1.0 - sumAlpha);
187
    for (phaseModel& phase : phases_)
188 189 190 191 192
    {
        volScalarField& alpha = phase;
        alpha += alpha*sumCorr;
    }

Henry's avatar
Henry committed
193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243
    calcAlphas();
}


Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseSystem::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::multiphaseSystem::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::multiphaseSystem::correctContactAngle
(
    const phaseModel& phase1,
    const phaseModel& phase2,
244
    surfaceVectorField::Boundary& nHatb
Henry's avatar
Henry committed
245 246
) const
{
247
    const volScalarField::Boundary& gbf
Henry's avatar
Henry committed
248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266
        = phase1.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]
            );

267 268
            const auto tp =
                acap.thetaProps().cfind(interfacePair(phase1, phase2));
Henry's avatar
Henry committed
269

270
            if (!tp.found())
Henry's avatar
Henry committed
271
            {
272 273
                FatalErrorInFunction
                    << "Cannot find interface " << interfacePair(phase1, phase2)
Henry's avatar
Henry committed
274 275 276 277 278 279 280
                    << "\n    in table of theta properties for patch "
                    << acap.patch().name()
                    << exit(FatalError);
            }

            bool matched = (tp.key().first() == phase1.name());

281
            const scalar theta0 = degToRad(tp().theta0(matched));
Henry's avatar
Henry committed
282 283 284 285 286 287 288
            scalarField theta(boundary[patchi].size(), theta0);

            scalar uTheta = tp().uTheta();

            // Calculate the dynamic contact angle if required
            if (uTheta > SMALL)
            {
289 290
                const scalar thetaA = degToRad(tp().thetaA(matched));
                const scalar thetaR = degToRad(tp().thetaR(matched));
Henry's avatar
Henry committed
291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350

                // Calculated the component of the velocity parallel to the wall
                vectorField Uwall
                (
                    phase1.U().boundaryField()[patchi].patchInternalField()
                  - phase1.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::multiphaseSystem::K
(
    const phaseModel& phase1,
    const phaseModel& phase2
) const
{
    tmp<surfaceVectorField> tnHatfv = nHatfv(phase1, phase2);

351
    correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryFieldRef());
Henry's avatar
Henry committed
352 353 354 355 356 357 358 359 360 361

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


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

Foam::multiphaseSystem::multiphaseSystem
(
Henry's avatar
Henry committed
362
    const volVectorField& U,
Henry's avatar
Henry committed
363 364 365
    const surfaceScalarField& phi
)
:
366 367 368 369 370 371 372 373 374 375 376 377
    IOdictionary
    (
        IOobject
        (
            "transportProperties",
            U.time().constant(),
            U.db(),
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
        )
    ),

Henry's avatar
Henry committed
378
    phases_(lookup("phases"), phaseModel::iNew(U.mesh())),
Henry's avatar
Henry committed
379

Henry's avatar
Henry committed
380
    mesh_(U.mesh()),
Henry's avatar
Henry committed
381 382 383 384 385 386 387 388 389 390 391 392 393
    phi_(phi),

    alphas_
    (
        IOobject
        (
            "alphas",
            mesh_.time().timeName(),
            mesh_,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh_,
394
        dimensionedScalar("alphas", dimless, 0.0)
Henry's avatar
Henry committed
395 396 397 398 399 400 401 402 403
    ),

    sigmas_(lookup("sigmas")),
    dimSigma_(1, 0, -2, 0, 0),
    cAlphas_(lookup("interfaceCompression")),
    Cvms_(lookup("virtualMass")),
    deltaN_
    (
        "deltaN",
404
        1e-8/cbrt(average(mesh_.V()))
Henry's avatar
Henry committed
405 406 407 408 409 410 411
    )
{
    calcAlphas();
    alphas_.write();

    interfaceDictTable dragModelsDict(lookup("drag"));

412
    forAllConstIters(dragModelsDict, iter)
Henry's avatar
Henry committed
413
    {
414
        dragModels_.set
Henry's avatar
Henry committed
415 416 417 418 419
        (
            iter.key(),
            dragModel::New
            (
                iter(),
420 421
                *phases_.lookup(iter.key().first()),
                *phases_.lookup(iter.key().second())
422
            ).ptr()
Henry's avatar
Henry committed
423 424
        );
    }
425

426
    for (const phaseModel& phase1 : phases_)
427
    {
428
        for (const phaseModel& phase2 : phases_)
429
        {
430
            if (&phase2 == &phase1)
431
            {
432 433
                continue;
            }
434

435 436 437 438 439 440 441 442 443
            const interfacePair key(phase1, phase2);

            if (sigmas_.found(key) && !cAlphas_.found(key))
            {
                WarningInFunction
                    << "Compression coefficient not specified for phase pair ("
                    << phase1.name() << ' ' << phase2.name()
                    << ") for which a surface tension coefficient is specified"
                    << endl;
444 445 446
            }
        }
    }
Henry's avatar
Henry committed
447 448 449 450 451 452 453
}


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

Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::rho() const
{
454
    auto iter = phases_.cbegin();
Henry's avatar
Henry committed
455 456

    tmp<volScalarField> trho = iter()*iter().rho();
457
    volScalarField& rho = trho.ref();
Henry's avatar
Henry committed
458

459
    for (++iter; iter != phases_.cend(); ++iter)
Henry's avatar
Henry committed
460
    {
461
        rho += iter()*iter().rho();
Henry's avatar
Henry committed
462 463 464 465 466 467
    }

    return trho;
}


468 469 470
Foam::tmp<Foam::scalarField>
Foam::multiphaseSystem::rho(const label patchi) const
{
471
    auto iter = phases_.cbegin();
472 473

    tmp<scalarField> trho = iter().boundaryField()[patchi]*iter().rho().value();
474
    scalarField& rho = trho.ref();
475

476
    for (++iter; iter != phases_.cend(); ++iter)
477
    {
478
        rho += iter().boundaryField()[patchi]*iter().rho().value();
479 480 481 482 483 484
    }

    return trho;
}


Henry's avatar
Henry committed
485 486
Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::nu() const
{
487
    auto iter = phases_.cbegin();
Henry's avatar
Henry committed
488

489
    tmp<volScalarField> tmu = iter()*(iter().rho()*iter().nu());
490
    volScalarField& mu = tmu.ref();
491

492
    for (++iter; iter != phases_.cend(); ++iter)
493
    {
494
        mu += iter()*(iter().rho()*iter().nu());
495 496 497 498 499 500 501 502 503
    }

    return tmu/rho();
}


Foam::tmp<Foam::scalarField>
Foam::multiphaseSystem::nu(const label patchi) const
{
504
    auto iter = phases_.cbegin();
505 506 507 508

    tmp<scalarField> tmu =
        iter().boundaryField()[patchi]
       *(iter().rho().value()*iter().nu().value());
509
    scalarField& mu = tmu.ref();
Henry's avatar
Henry committed
510

511
    for (++iter; iter != phases_.cend(); ++iter)
Henry's avatar
Henry committed
512
    {
513
        mu +=
514 515
            iter().boundaryField()[patchi]
           *(iter().rho().value()*iter().nu().value());
Henry's avatar
Henry committed
516 517
    }

518
    return tmu/rho(patchi);
Henry's avatar
Henry committed
519 520 521
}


Henry's avatar
Henry committed
522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537
Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::Cvm
(
    const phaseModel& phase
) const
{
    tmp<volScalarField> tCvm
    (
        new volScalarField
        (
            IOobject
            (
                "Cvm",
                mesh_.time().timeName(),
                mesh_
            ),
            mesh_,
538
            dimensionedScalar(dimDensity, Zero)
Henry's avatar
Henry committed
539 540 541
        )
    );

542
    for (const phaseModel& phase2 : phases_)
Henry's avatar
Henry committed
543
    {
544 545 546 547 548 549
        if (&phase2 == &phase)
        {
            continue;
        }

        auto iterCvm = Cvms_.cfind(interfacePair(phase, phase2));
Henry's avatar
Henry committed
550

551
        if (iterCvm.found())
Henry's avatar
Henry committed
552
        {
553 554 555 556 557
            tCvm.ref() += iterCvm()*phase2.rho()*phase2;
        }
        else
        {
            iterCvm = Cvms_.cfind(interfacePair(phase2, phase));
558

559
            if (iterCvm.found())
560
            {
561
                tCvm.ref() += iterCvm()*phase.rho()*phase2;
562
            }
Henry's avatar
Henry committed
563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585
        }
    }

    return tCvm;
}


Foam::tmp<Foam::volVectorField> Foam::multiphaseSystem::Svm
(
    const phaseModel& phase
) const
{
    tmp<volVectorField> tSvm
    (
        new volVectorField
        (
            IOobject
            (
                "Svm",
                mesh_.time().timeName(),
                mesh_
            ),
            mesh_,
586 587 588 589 590 591
            dimensionedVector
            (
                "Svm",
                dimensionSet(1, -2, -2, 0, 0),
                Zero
            )
Henry's avatar
Henry committed
592 593 594
        )
    );

595
    for (const phaseModel& phase2 : phases_)
Henry's avatar
Henry committed
596
    {
597 598 599 600 601 602
        if (&phase2 == &phase)
        {
            continue;
        }

        auto Cvm = Cvms_.cfind(interfacePair(phase, phase2));
Henry's avatar
Henry committed
603

604
        if (Cvm.found())
Henry's avatar
Henry committed
605
        {
606 607 608 609 610
            tSvm.ref() += Cvm()*phase2.rho()*phase2*phase2.DDtU();
        }
        else
        {
            Cvm = Cvms_.cfind(interfacePair(phase2, phase));
611

612
            if (Cvm.found())
613
            {
614
                tSvm.ref() += Cvm()*phase.rho()*phase2*phase2.DDtU();
615
            }
Henry's avatar
Henry committed
616 617 618
        }
    }

619
    volVectorField::Boundary& SvmBf =
620 621
        tSvm.ref().boundaryFieldRef();

622
    // Remove virtual mass at fixed-flux boundaries
623 624 625 626 627 628 629 630 631 632
    forAll(phase.phi().boundaryField(), patchi)
    {
        if
        (
            isA<fixedValueFvsPatchScalarField>
            (
                phase.phi().boundaryField()[patchi]
            )
        )
        {
633
            SvmBf[patchi] = Zero;
634 635 636
        }
    }

Henry's avatar
Henry committed
637 638 639 640 641 642 643
    return tSvm;
}


Foam::autoPtr<Foam::multiphaseSystem::dragCoeffFields>
Foam::multiphaseSystem::dragCoeffs() const
{
644
    autoPtr<dragCoeffFields> dragCoeffsPtr(new dragCoeffFields);
Henry's avatar
Henry committed
645

646
    forAllConstIters(dragModels_, iter)
Henry's avatar
Henry committed
647 648 649
    {
        const dragModel& dm = *iter();

650
        volScalarField* Kptr =
Henry's avatar
Henry committed
651
            (
652 653
                max
                (
654 655
                    // fvc::average(dm.phase1()*dm.phase2()),
                    // fvc::average(dm.phase1())*fvc::average(dm.phase2()),
656
                    dm.phase1()*dm.phase2(),
657 658 659 660 661 662 663 664 665 666
                    dm.residualPhaseFraction()
                )
               *dm.K
                (
                    max
                    (
                        mag(dm.phase1().U() - dm.phase2().U()),
                        dm.residualSlip()
                    )
                )
667 668
            ).ptr();

669
        volScalarField::Boundary& Kbf = Kptr->boundaryFieldRef();
670

671 672 673 674 675 676 677 678 679 680 681
        // Remove drag at fixed-flux boundaries
        forAll(dm.phase1().phi().boundaryField(), patchi)
        {
            if
            (
                isA<fixedValueFvsPatchScalarField>
                (
                    dm.phase1().phi().boundaryField()[patchi]
                )
            )
            {
682
                Kbf[patchi] = 0.0;
683 684 685
            }
        }

686
        dragCoeffsPtr().set(iter.key(), Kptr);
Henry's avatar
Henry committed
687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709
    }

    return dragCoeffsPtr;
}


Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::dragCoeff
(
    const phaseModel& phase,
    const dragCoeffFields& dragCoeffs
) const
{
    tmp<volScalarField> tdragCoeff
    (
        new volScalarField
        (
            IOobject
            (
                "dragCoeff",
                mesh_.time().timeName(),
                mesh_
            ),
            mesh_,
710 711 712 713 714 715
            dimensionedScalar
            (
                "dragCoeff",
                dimensionSet(1, -3, -1, 0, 0),
                0
            )
Henry's avatar
Henry committed
716 717 718 719 720 721 722 723
        )
    );

    dragModelTable::const_iterator dmIter = dragModels_.begin();
    dragCoeffFields::const_iterator dcIter = dragCoeffs.begin();
    for
    (
        ;
724
        dmIter.good() && dcIter.good();
Henry's avatar
Henry committed
725 726 727 728 729 730 731 732 733
        ++dmIter, ++dcIter
    )
    {
        if
        (
            &phase == &dmIter()->phase1()
         || &phase == &dmIter()->phase2()
        )
        {
734
            tdragCoeff.ref() += *dcIter();
Henry's avatar
Henry committed
735 736 737 738 739 740 741
        }
    }

    return tdragCoeff;
}


742 743 744 745
Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseSystem::surfaceTension
(
    const phaseModel& phase1
) const
Henry's avatar
Henry committed
746
{
747
    tmp<surfaceScalarField> tSurfaceTension
Henry's avatar
Henry committed
748 749 750 751 752
    (
        new surfaceScalarField
        (
            IOobject
            (
753
                "surfaceTension",
Henry's avatar
Henry committed
754 755 756 757
                mesh_.time().timeName(),
                mesh_
            ),
            mesh_,
758 759 760 761 762 763
            dimensionedScalar
            (
                "surfaceTension",
                dimensionSet(1, -2, -2, 0, 0),
                0
            )
Henry's avatar
Henry committed
764 765
        )
    );
766
    tSurfaceTension.ref().setOriented();
767

768
    for (const phaseModel& phase2 : phases_)
Henry's avatar
Henry committed
769
    {
770
        if (&phase2 == &phase1)
Henry's avatar
Henry committed
771
        {
772 773
            continue;
        }
Henry's avatar
Henry committed
774

775 776 777 778 779 780 781 782 783 784 785
        const auto sigma = sigmas_.cfind(interfacePair(phase1, phase2));

        if (sigma.found())
        {
            tSurfaceTension.ref() +=
                dimensionedScalar("sigma", dimSigma_, *sigma)
               *fvc::interpolate(K(phase1, phase2))*
                (
                    fvc::interpolate(phase2)*fvc::snGrad(phase1)
                  - fvc::interpolate(phase1)*fvc::snGrad(phase2)
                );
Henry's avatar
Henry committed
786 787 788
        }
    }

789
    return tSurfaceTension;
Henry's avatar
Henry committed
790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806
}


Foam::tmp<Foam::volScalarField>
Foam::multiphaseSystem::nearInterface() const
{
    tmp<volScalarField> tnearInt
    (
        new volScalarField
        (
            IOobject
            (
                "nearInterface",
                mesh_.time().timeName(),
                mesh_
            ),
            mesh_,
807
            dimensionedScalar("nearInterface", dimless, 0.0)
Henry's avatar
Henry committed
808 809 810
        )
    );

811
    for (const phaseModel& phase : phases_)
Henry's avatar
Henry committed
812
    {
Henry Weller's avatar
Henry Weller committed
813
        tnearInt.ref() =
814
            max(tnearInt(), pos0(phase - 0.01)*pos0(0.99 - phase));
Henry's avatar
Henry committed
815 816 817 818 819 820 821 822
    }

    return tnearInt;
}


void Foam::multiphaseSystem::solve()
{
823
    for (phaseModel& phase : phases_)
Henry's avatar
Henry committed
824
    {
825
        phase.correct();
Henry's avatar
Henry committed
826 827 828 829
    }

    const Time& runTime = mesh_.time();

830
    const dictionary& alphaControls = mesh_.solverDict("alpha");
831
    label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
Henry's avatar
Henry committed
832 833 834 835 836

    if (nAlphaSubCycles > 1)
    {
        dimensionedScalar totalDeltaT = runTime.deltaT();

Henry's avatar
Henry committed
837
        PtrList<volScalarField> alpha0s(phases_.size());
838
        PtrList<surfaceScalarField> alphaPhiSums(phases_.size());
Henry's avatar
Henry committed
839

840 841
        label phasei = 0;
        for (phaseModel& phase : phases_)
Henry's avatar
Henry committed
842 843 844 845 846 847 848 849 850
        {
            volScalarField& alpha = phase;

            alpha0s.set
            (
                phasei,
                new volScalarField(alpha.oldTime())
            );

851
            alphaPhiSums.set
Henry's avatar
Henry committed
852 853 854 855 856 857 858 859 860 861 862
            (
                phasei,
                new surfaceScalarField
                (
                    IOobject
                    (
                        "phiSum" + alpha.name(),
                        runTime.timeName(),
                        mesh_
                    ),
                    mesh_,
863
                    dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0)
Henry's avatar
Henry committed
864 865 866
                )
            );

867
            ++phasei;
Henry's avatar
Henry committed
868 869
        }

Henry's avatar
Henry committed
870 871
        for
        (
Henry's avatar
Henry committed
872 873 874 875 876
            subCycleTime alphaSubCycle
            (
                const_cast<Time&>(runTime),
                nAlphaSubCycles
            );
Henry's avatar
Henry committed
877 878 879 880
            !(++alphaSubCycle).end();
        )
        {
            solveAlphas();
Henry's avatar
Henry committed
881

882 883
            label phasei = 0;
            for (const phaseModel& phase : phases_)
Henry's avatar
Henry committed
884
            {
885 886 887
                alphaPhiSums[phasei] += phase.alphaPhi()/nAlphaSubCycles;

                ++phasei;
Henry's avatar
Henry committed
888 889 890 891
            }
        }

        phasei = 0;
892
        for (phaseModel& phase : phases_)
Henry's avatar
Henry committed
893 894 895
        {
            volScalarField& alpha = phase;

896
            phase.alphaPhi() = alphaPhiSums[phasei];
Henry's avatar
Henry committed
897

Henry's avatar
Henry committed
898 899
            // Correct the time index of the field
            // to correspond to the global time
Henry's avatar
Henry committed
900 901 902 903 904 905
            alpha.timeIndex() = runTime.timeIndex();

            // Reset the old-time field value
            alpha.oldTime() = alpha0s[phasei];
            alpha.oldTime().timeIndex() = runTime.timeIndex();

906
            ++phasei;
Henry's avatar
Henry committed
907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924
        }
    }
    else
    {
        solveAlphas();
    }
}


bool Foam::multiphaseSystem::read()
{
    if (regIOobject::read())
    {
        bool readOK = true;

        PtrList<entry> phaseData(lookup("phases"));
        label phasei = 0;

925
        for (phaseModel& phase : phases_)
Henry's avatar
Henry committed
926
        {
927
            readOK &= phase.read(phaseData[phasei++].dict());
Henry's avatar
Henry committed
928 929
        }

930 931 932
        lookup("sigmas") >> sigmas_;
        lookup("interfaceCompression") >> cAlphas_;
        lookup("virtualMass") >> Cvms_;
Henry's avatar
Henry committed
933 934 935

        return readOK;
    }
936 937 938 939
    else
    {
        return false;
    }
Henry's avatar
Henry committed
940 941 942 943
}


// ************************************************************************* //