multiphaseSystem.C 18.6 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-2019 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
-------------------------------------------------------------------------------
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"

#include "MULES.H"
#include "subCycle.H"
34
#include "UniformField.H"
35 36 37 38 39 40 41 42 43 44 45 46

#include "fvcDdt.H"
#include "fvcDiv.H"
#include "fvcSnGrad.H"
#include "fvcFlux.H"
#include "fvcMeshPhi.H"
#include "fvcSup.H"

#include "fvmDdt.H"
#include "fvmLaplacian.H"
#include "fvmSup.H"

47 48
#include "unitConversion.H"

49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
// * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //

namespace Foam
{
    defineTypeNameAndDebug(multiphaseSystem, 0);
    defineRunTimeSelectionTable(multiphaseSystem, dictionary);
}


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

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

sergio's avatar
sergio committed
65
    forAll(phases(), i)
66
    {
sergio's avatar
sergio committed
67
        alphas_ += level*phases()[i];
68 69 70 71 72 73 74
        level += 1.0;
    }
}


void Foam::multiphaseSystem::solveAlphas()
{
sergio's avatar
sergio committed
75
    forAll(phases(), phasei)
76
    {
sergio's avatar
sergio committed
77
        phases()[phasei].correctBoundaryConditions();
78 79
    }

80 81 82 83 84 85 86 87 88 89 90 91
    // Calculate the void fraction
    volScalarField alphaVoid
    (
        IOobject
        (
            "alphaVoid",
            mesh_.time().timeName(),
            mesh_
        ),
        mesh_,
        dimensionedScalar("one", dimless, 1)
    );
sergio's avatar
sergio committed
92
    forAll(stationaryPhases(), stationaryPhasei)
93
    {
sergio's avatar
sergio committed
94
        alphaVoid -= stationaryPhases()[stationaryPhasei];
95
    }
96

97 98 99
    // Generate face-alphas
    PtrList<surfaceScalarField> alphafs(phases().size());
    forAll(phases(), phasei)
100
    {
101 102 103 104 105 106 107 108 109 110 111
        phaseModel& phase = phases()[phasei];
        alphafs.set
        (
            phasei,
            new surfaceScalarField
            (
                IOobject::groupName("alphaf", phase.name()),
                upwind<scalar>(mesh_, phi_).interpolate(phase)
            )
        );
    }
112

113 114
    // Create correction fluxes
    PtrList<surfaceScalarField> alphaPhiCorrs(phases().size());
sergio's avatar
sergio committed
115
    forAll(stationaryPhases(), stationaryPhasei)
116
    {
sergio's avatar
sergio committed
117 118
        phaseModel& phase = stationaryPhases()[stationaryPhasei];

119 120
        alphaPhiCorrs.set
        (
121
            phase.index(),
122 123
            new surfaceScalarField
            (
124 125 126 127 128
                IOobject::groupName("alphaPhiCorr", phase.name()),
              - upwind<scalar>(mesh_, phi_).flux(phase)
            )
        );
    }
sergio's avatar
sergio committed
129
    forAll(movingPhases(), movingPhasei)
130
    {
sergio's avatar
sergio committed
131 132
        phaseModel& phase = movingPhases()[movingPhasei];
        volScalarField& alpha = phase;
133 134 135 136 137 138 139 140

        alphaPhiCorrs.set
        (
            phase.index(),
            new surfaceScalarField
            (
                IOobject::groupName("alphaPhiCorr", phase.name()),
                fvc::flux(phi_, phase, "div(phi," + alpha.name() + ')')
141 142 143
            )
        );

144
        surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phase.index()];
145

sergio's avatar
sergio committed
146
        forAll(phases(), phasei)
147
        {
sergio's avatar
sergio committed
148 149
            phaseModel& phase2 = phases()[phasei];
            volScalarField& alpha2 = phase2;
150 151 152 153 154

            if (&phase2 == &phase) continue;

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

sergio's avatar
sergio committed
155 156 157 158
            cAlphaTable::const_iterator cAlpha
            (
                cAlphas_.find(phasePairKey(phase.name(), phase2.name()))
            );
159

sergio's avatar
sergio committed
160
            if (cAlpha != cAlphas_.end())
161 162 163 164 165 166 167 168 169
            {
                surfaceScalarField phic
                (
                    (mag(phi_) + mag(phir))/mesh_.magSf()
                );

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

sergio's avatar
sergio committed
170
            word phirScheme
171
            (
172
                "div(phir," + alpha2.name() + ',' + alpha.name() + ')'
173 174 175 176 177 178 179 180 181 182
            );

            alphaPhiCorr += fvc::flux
            (
                -fvc::flux(-phir, phase2, phirScheme),
                phase,
                phirScheme
            );
        }

183
        phase.correctInflowOutflow(alphaPhiCorr);
184

185
        MULES::limit
186
        (
187 188 189 190 191 192 193 194 195 196 197
            geometricOneField(),
            phase,
            phi_,
            alphaPhiCorr,
            zeroField(),
            zeroField(),
            min(alphaVoid.primitiveField(), phase.alphaMax())(),
            zeroField(),
            true
        );
    }
198

199 200
    // Limit the flux sums, fixing those of the stationary phases
    labelHashSet fixedAlphaPhiCorrs;
sergio's avatar
sergio committed
201
    forAll(stationaryPhases(), stationaryPhasei)
202
    {
sergio's avatar
sergio committed
203
        fixedAlphaPhiCorrs.insert(stationaryPhases()[stationaryPhasei].index());
204 205
    }
    MULES::limitSum(alphafs, alphaPhiCorrs, fixedAlphaPhiCorrs);
206

207
    // Solve for the moving phase alphas
sergio's avatar
sergio committed
208
    forAll(movingPhases(), movingPhasei)
209
    {
sergio's avatar
sergio committed
210
        phaseModel& phase = movingPhases()[movingPhasei];
211 212
        volScalarField& alpha = phase;

213
        surfaceScalarField& alphaPhi = alphaPhiCorrs[phase.index()];
214 215
        alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase);
        phase.correctInflowOutflow(alphaPhi);
216

217
        volScalarField::Internal Sp
218 219 220 221 222 223 224 225
        (
            IOobject
            (
                "Sp",
                mesh_.time().timeName(),
                mesh_
            ),
            mesh_,
sergio's avatar
sergio committed
226
            dimensionedScalar("zero", dimless/dimTime, 0)
227 228
        );

229
        volScalarField::Internal Su
230
        (
231 232
            "Su",
            min(alpha, scalar(1))*fvc::div(fvc::absolute(phi_, phase.U()))
233 234
        );

235
        if (phase.divU().valid())
236
        {
237
            const scalarField& dgdt = phase.divU()();
238 239 240 241 242 243 244 245 246 247 248 249

            forAll(dgdt, celli)
            {
                if (dgdt[celli] > 0.0)
                {
                    Sp[celli] -= dgdt[celli];
                    Su[celli] += dgdt[celli];
                }
                else if (dgdt[celli] < 0.0)
                {
                    Sp[celli] +=
                        dgdt[celli]
250
                       *(1 - alpha[celli])/max(alpha[celli], 1e-4);
251 252 253 254
                }
            }
        }

sergio's avatar
sergio committed
255
        forAll(phases(), phasej)
256
        {
sergio's avatar
sergio committed
257
            const phaseModel& phase2 = phases()[phasej];
258 259 260 261
            const volScalarField& alpha2 = phase2;

            if (&phase2 == &phase) continue;

262
            if (phase2.divU().valid())
263
            {
264
                const scalarField& dgdt2 = phase2.divU()();
265 266 267 268 269 270 271

                forAll(dgdt2, celli)
                {
                    if (dgdt2[celli] < 0.0)
                    {
                        Sp[celli] +=
                            dgdt2[celli]
272
                           *(1 - alpha2[celli])/max(alpha2[celli], 1e-4);
273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289

                        Su[celli] -=
                            dgdt2[celli]
                           *alpha[celli]/max(alpha2[celli], 1e-4);
                    }
                    else if (dgdt2[celli] > 0.0)
                    {
                        Sp[celli] -= dgdt2[celli];
                    }
                }
            }
        }

        MULES::explicitSolve
        (
            geometricOneField(),
            alpha,
290
            alphaPhi,
291 292 293 294
            Sp,
            Su
        );

295 296
        phase.alphaPhiRef() = alphaPhi;
    }
297

298
    // Report the phase fractions and the phase fraction sum
sergio's avatar
sergio committed
299
    forAll(phases(), phasei)
300
    {
sergio's avatar
sergio committed
301 302
        phaseModel& phase = phases()[phasei];

303
        Info<< phase.name() << " fraction, min, max = "
304 305 306 307
            << phase.weightedAverage(mesh_.V()).value()
            << ' ' << min(phase).value()
            << ' ' << max(phase).value()
            << endl;
308
    }
309

310 311 312 313 314 315 316 317 318
    volScalarField sumAlphaMoving
    (
        IOobject
        (
            "sumAlphaMoving",
            mesh_.time().timeName(),
            mesh_
        ),
        mesh_,
sergio's avatar
sergio committed
319
        dimensionedScalar("zero", dimless, 0)
320
    );
sergio's avatar
sergio committed
321
    forAll(movingPhases(), movingPhasei)
322
    {
sergio's avatar
sergio committed
323
        sumAlphaMoving += movingPhases()[movingPhasei];
324 325 326
    }

    Info<< "Phase-sum volume fraction, min, max = "
327 328 329
        << (sumAlphaMoving + 1 - alphaVoid)().weightedAverage(mesh_.V()).value()
        << ' ' << min(sumAlphaMoving + 1 - alphaVoid).value()
        << ' ' << max(sumAlphaMoving + 1 - alphaVoid).value()
330
        << endl;
331

332
    // Correct the sum of the phase fractions to avoid drift
sergio's avatar
sergio committed
333
    forAll(movingPhases(), movingPhasei)
334
    {
sergio's avatar
sergio committed
335
        movingPhases()[movingPhasei] *= alphaVoid/sumAlphaMoving;
336
    }
337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386
}


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,
387
    surfaceVectorField::Boundary& nHatb
388 389
) const
{
390
    const volScalarField::Boundary& gbf
391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409
        = 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]
            );

sergio's avatar
sergio committed
410 411 412 413
            alphaContactAngleFvPatchScalarField::thetaPropsTable::
                const_iterator tp =
                    acap.thetaProps()
                   .find(phasePairKey(phase1.name(), phase2.name()));
414

sergio's avatar
sergio committed
415
            if (tp == acap.thetaProps().end())
416
            {
417 418
                FatalErrorInFunction
                    << "Cannot find interface "
419 420 421 422 423 424
                    << phasePairKey(phase1.name(), phase2.name())
                    << "\n    in table of theta properties for patch "
                    << acap.patch().name()
                    << exit(FatalError);
            }

sergio's avatar
sergio committed
425
            bool matched = (tp.key().first() == phase1.name());
426

427
            scalar theta0 = degToRad(tp().theta0(matched));
428 429 430 431 432 433 434
            scalarField theta(boundary[patchi].size(), theta0);

            scalar uTheta = tp().uTheta();

            // Calculate the dynamic contact angle if required
            if (uTheta > SMALL)
            {
435 436
                const scalar thetaA = degToRad(tp().thetaA(matched));
                const scalar thetaR = degToRad(tp().thetaR(matched));
437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475

                // 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]);
            }

476
            scalarField det(1 - a12*a12);
477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496

            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);

497
    correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryFieldRef());
498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517

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


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

Foam::multiphaseSystem::multiphaseSystem
(
    const fvMesh& mesh
)
:
    phaseSystem(mesh),

    alphas_
    (
        IOobject
        (
            "alphas",
518
            mesh_.time().timeName(),
519 520 521 522 523
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh,
sergio's avatar
sergio committed
524
        dimensionedScalar("zero", dimless, 0)
525 526 527 528 529 530 531
    ),

    cAlphas_(lookup("interfaceCompression")),

    deltaN_
    (
        "deltaN",
sergio's avatar
sergio committed
532
        1e-8/pow(average(mesh_.V()), 1.0/3.0)
533 534
    )
{
sergio's avatar
sergio committed
535
    forAll(phases(), phasei)
536
    {
sergio's avatar
sergio committed
537 538
        volScalarField& alphai = phases()[phasei];
        mesh_.setFluxRequired(alphai.name());
539 540 541 542
    }
}


sergio's avatar
sergio committed
543 544 545 546 547 548
// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::multiphaseSystem::~multiphaseSystem()
{}


549 550 551 552 553 554 555 556 557
// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //

Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseSystem::surfaceTension
(
    const phaseModel& phase1
) const
{
    tmp<surfaceScalarField> tSurfaceTension
    (
558
        surfaceScalarField::New
559
        (
560
            "surfaceTension",
561
            mesh_,
sergio's avatar
sergio committed
562
            dimensionedScalar("zero", dimensionSet(1, -2, -2, 0, 0), 0)
563 564 565
        )
    );

566
    tSurfaceTension.ref().setOriented();
567

sergio's avatar
sergio committed
568
    forAll(phases(), phasej)
569
    {
sergio's avatar
sergio committed
570 571
        const phaseModel& phase2 = phases()[phasej];

572 573 574 575
        if (&phase2 != &phase1)
        {
            phasePairKey key12(phase1.name(), phase2.name());

sergio's avatar
sergio committed
576
            cAlphaTable::const_iterator cAlpha(cAlphas_.find(key12));
577

sergio's avatar
sergio committed
578
            if (cAlpha != cAlphas_.end())
579
            {
580
                tSurfaceTension.ref() +=
581 582 583 584 585 586 587 588
                    fvc::interpolate(sigma(key12)*K(phase1, phase2))
                   *(
                        fvc::interpolate(phase2)*fvc::snGrad(phase1)
                      - fvc::interpolate(phase1)*fvc::snGrad(phase2)
                    );
            }
        }
    }
sergio's avatar
sergio committed
589

590
    tSurfaceTension->setOriented();
591 592 593 594 595 596 597 598 599 600

    return tSurfaceTension;
}


Foam::tmp<Foam::volScalarField>
Foam::multiphaseSystem::nearInterface() const
{
    tmp<volScalarField> tnearInt
    (
601
        volScalarField::New
602
        (
603
            "nearInterface",
604
            mesh_,
sergio's avatar
sergio committed
605
            dimensionedScalar("zero", dimless, 0)
606 607 608
        )
    );

sergio's avatar
sergio committed
609
    forAll(phases(), phasei)
610
    {
611
        tnearInt.ref() = max
612 613
        (
            tnearInt(),
sergio's avatar
sergio committed
614
            pos0(phases()[phasei] - 0.01)*pos0(0.99 - phases()[phasei])
615
        );
616 617 618 619 620 621 622 623
    }

    return tnearInt;
}


void Foam::multiphaseSystem::solve()
{
624
    const Time& runTime = mesh_.time();
625 626

    const dictionary& alphaControls = mesh_.solverDict("alpha");
627
    label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
628

629
    bool LTS = fv::localEulerDdt::enabled(mesh_);
630 631 632 633 634 635 636 637

    if (nAlphaSubCycles > 1)
    {
        tmp<volScalarField> trSubDeltaT;

        if (LTS)
        {
            trSubDeltaT =
638
                fv::localEulerDdt::localRSubDeltaT(mesh_, nAlphaSubCycles);
639 640
        }

641
        List<volScalarField*> alphaPtrs(phases().size());
642
        PtrList<surfaceScalarField> alphaPhiSums(phases().size());
643

sergio's avatar
sergio committed
644
        forAll(phases(), phasei)
645
        {
sergio's avatar
sergio committed
646
            phaseModel& phase = phases()[phasei];
647 648
            volScalarField& alpha = phase;

649
            alphaPtrs[phasei] = &alpha;
650

651
            alphaPhiSums.set
652 653 654 655 656 657 658 659 660 661 662
            (
                phasei,
                new surfaceScalarField
                (
                    IOobject
                    (
                        "phiSum" + alpha.name(),
                        runTime.timeName(),
                        mesh_
                    ),
                    mesh_,
sergio's avatar
sergio committed
663
                    dimensionedScalar("zero", dimensionSet(0, 3, -1, 0, 0), 0)
664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679
                )
            );
        }

        for
        (
            subCycleTime alphaSubCycle
            (
                const_cast<Time&>(runTime),
                nAlphaSubCycles
            );
            !(++alphaSubCycle).end();
        )
        {
            solveAlphas();

sergio's avatar
sergio committed
680
            forAll(phases(), phasei)
681
            {
sergio's avatar
sergio committed
682
                alphaPhiSums[phasei] += phases()[phasei].alphaPhi();
683 684 685
            }
        }

sergio's avatar
sergio committed
686
        forAll(phases(), phasei)
687
        {
sergio's avatar
sergio committed
688
            phaseModel& phase = phases()[phasei];
689
            if (phase.stationary()) continue;
690

691
            phase.alphaPhiRef() = alphaPhiSums[phasei]/nAlphaSubCycles;
692 693 694 695 696 697 698
        }
    }
    else
    {
        solveAlphas();
    }

sergio's avatar
sergio committed
699
    forAll(phases(), phasei)
700
    {
sergio's avatar
sergio committed
701
        phaseModel& phase = phases()[phasei];
702 703 704 705
        if (phase.stationary()) continue;

        phase.alphaRhoPhiRef() =
            fvc::interpolate(phase.rho())*phase.alphaPhi();
706

707
        phase.clip(SMALL, 1 - SMALL);
708 709 710 711 712 713 714
    }

    calcAlphas();
}


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