[tgvosk]

Hello,

I have a question for any of the MULES experts out there.

I am attempting to add source terms to multiphaseInterFoam and am running into trouble maintaining boundedness. I have added the source term for each phase to both the limiter for that phase (MULES::limit) and the solver for that phase (MULES::explicitSolve) but the solution still slowly becomes unbounded when the source terms are non-negligible. Odder still, increasing nAlphaSubCycles only makes it become unbounded sooner.

Does anyone know what could cause MULES to stop producing bounded results? It seems like the limiter algorithm called on phase A isn't being influenced enough by the source in phase B, for example, but I'm not sure how that could be fixed.

Thanks!

[hrushi.397]

Hi Tyler,

Did you find solution to your problem?

Thanks and regards,

Hrushi

[tgvosk]

Hrushi,

I never did find a solution to this, so I ended up restructuring my code to use only two phases (liquid and gas) and differentiate between components using species mass fraction instead. With only two phases it's easier to maintain boundedness and the MULES multiphase limiter gets really expensive with more phases anyway.

-Tyler

[hrushi.397]

Thank you for a quick response, Tyler.

I am facing a similar issue with two phase simulations. Is it possible to post your implementation of source for two phases for me to compare?

Thanks and regards,

Hrushi

[tgvosk]

With only two phases, there are several solvers that do this already. Take a look at the alphaEqn.H files in compressibleInterFoam and interPhaseChangeFoam for examples. If you have severe boundedness issues, the implicit solver used in interPhaseChangeFoam may help.

-Tyler

[hrushi.397]

Thanks Tyler. Will look them up.

- Hrushi

[Villo]

Hi foamers!
Tyler, which OpenFoam release are you using?
I had similar problems for multiphase flow (twoLiquidMixingFoam and interMixingFoam) in OpenFoam 2.0. I have try hundreds different numerical schemes to limit the field, but problems were finally suppressed moving from OpenFoam 2.0 to OpenFoam 2.2.0 (more limiters (as lambda) and better structured)).
Regards

[Mastra]

Dear all,

I am also trying to implement a source term into the phaseEquation of a multiphase solver (multiphaseInterFoam).

My goal is to convert one phase e.g. alpha2 ( 3 phases in total: alpha1, alpha2, alpha3) into another phase (e.g. alpha3 with different transportProperties -> liquid into vapour) at certain conditions.
The distinction criterion for each phase should be a volScalarField between 0 and 1, where 0 means alpha2 only and 1 alpha3 only respectively.

Is this possible using MULES in the solveAlphas() Method in multiphaseMixture.C ?

void Foam::multiphaseMixture::solveAlphas
(
    const scalar cAlpha
)
{
    static label nSolves=-1;
    nSolves++;

    word alphaScheme("div(phi,alpha)");
    word alpharScheme("div(phirb,alpha)");

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

    PtrList<surfaceScalarField> phiAlphaCorrs(phases_.size());
    int phasei = 0;

    forAllIter(PtrDictionary<phase>, phases_, iter)
    {
        phase& alpha = iter();

        phiAlphaCorrs.set
        (
            phasei,
            new surfaceScalarField
            (
                fvc::flux
                (
                    phi_,
                    alpha,
                    alphaScheme
                )
            )
        );

        surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei];

        forAllIter(PtrDictionary<phase>, phases_, iter2)
        {
            phase& alpha2 = iter2();

            if (&alpha2 == &alpha) continue;

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

            phiAlphaCorr += fvc::flux
            (
                -fvc::flux(-phir, alpha2, alpharScheme),
                alpha,
                alpharScheme
            );
        }

        MULES::limit
        (
            geometricOneField(),
            alpha,
            phi_,
            phiAlphaCorr,
            zeroField(),
            zeroField(),
            1,
            0,
            3,
            true
        );

        phasei++;
    }

    MULES::limitSum(phiAlphaCorrs);

    rhoPhi_ = dimensionedScalar("0", dimensionSet(1, 0, -1, 0, 0), 0);

    volScalarField sumAlpha
    (
        IOobject
        (
            "sumAlpha",
            mesh_.time().timeName(),
            mesh_
        ),
        mesh_,
        dimensionedScalar("sumAlpha", dimless, 0)
    );

    phasei = 0;

    forAllIter(PtrDictionary<phase>, phases_, iter)
    {
        phase& alpha = iter();

        surfaceScalarField& phiAlpha = phiAlphaCorrs[phasei];
        phiAlpha += upwind<scalar>(mesh_, phi_).flux(alpha);

        MULES::explicitSolve
        (
            geometricOneField(),
            alpha,
            phiAlpha,
            zeroField(),
            zeroField()
        );

        rhoPhi_ += phiAlpha* fvc::interpolate (alpha.rho());

        Info<< alpha.name() << " volume fraction, min, max = "
            << alpha.weightedAverage(mesh_.V()).value()
            << ' ' << min(alpha).value()
            << ' ' << max(alpha).value()
            << endl;

        sumAlpha += alpha;

        phasei++;
    }

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

    calcAlphas();

}

I know i have to replace the zeroFields somehow with the implicit and explicit source terms, but how can i make sure that alpha2 converts into alpha3 only and alpha1 stays unchanged using the distinction Field between 0 and 1 ?

Or is there another way to do this ?

Thanks and Regards

Markus

[wavefunction]

Originally Posted by tgvosk

Hrushi,

I never did find a solution to this, so I ended up restructuring my code to use only two phases (liquid and gas) and differentiate between components using species mass fraction instead. With only two phases it's easier to maintain boundedness and the MULES multiphase limiter gets really expensive with more phases anyway.

-Tyler

Hi Tyler.

So if you modified interPhaseChangeFoam to accomodate your source terms in it, did you add an equation to solve alpha2 as well? Because if there are multiple phases present, alpha2 is not equal to (1 - alpha1) as given in the alphaEqn.H file in interPhaseChangeFoam.