[tom_flint2012]

Good afternoon,

I've been trying to add diffusion between selected phase pairs in multiphaseeulerfoam, in order to capture diffusion of phases during mixing, as the processes I am looking at are quite dependent on this.

I had some success in editing the solveAlphas() function to take a laplacian term as a source, but I'm afraid I dont know enough about the MULES algorithm to go much further, I could add a diffusive term to all phases but not to selected phases, at their interfaces.

I have created the constructs to read the diffusion coefficients for selected pairs (I call this dAlpha in the code) but when I try and apply the diffusion in this way there seems to be no effect.

The closest I have come to a solution is to create a volscalarfield of the diffusion constants, then calculate the diffusive flux based on this, however, in this way I get diffusion, not just at the desired pairs, and even if the selected pairs are not in contact.

Any help would be greatly appreciated. I tried to upload the solver but the forum complained that it was too large. If anyone would like to help with this I can send them the solver and hopefully we can crack this together!

void Foam::multiphaseSystem::solveAlphas() //const
{
    PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
    PtrList<volScalarField> alphalaps(phases_.size());




////////////////////////////////Tried to make some fields to keep diffusion in desired phases
volScalarField difffield
    (
        IOobject
        (
            "difffield",
            mesh_.time().timeName(),
            mesh_
        ),
        mesh_,
        dimensionedScalar("difffield", dimensionSet(0, 0, -1, 0, 0), 0)
    );


    volScalarField coefffield
    (
        IOobject
        (
            "coefffield",
            mesh_.time().timeName(),
            mesh_
        ),
        mesh_,
        dimensionedScalar("coefffield", dimensionSet(0, 2, -1, 0, 0), 0)
    );

    forAllIter(PtrDictionary<phaseModel>, phases_, iter)
    {
        phaseModel& phase1 = iter();
        volScalarField& alpha1 = phase1;

        coefffield+=alpha1*phase1.diffusionconstant();

        }
//////////////////////////////////////////////////////////////////////////////////////////


    int phasei = 0;

    forAllIter(PtrDictionary<phaseModel>, phases_, iter)
    {
        phaseModel& phase1 = iter();
        volScalarField& alpha1 = phase1;

        phase1.alphaPhi() =
            dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0);

        alphaPhiCorrs.set
        (
            phasei,
            new surfaceScalarField
            (
                "phi" + alpha1.name() + "Corr",
                fvc::flux
                (
                    phi_,
                    phase1,
                    "div(phi," + alpha1.name() + ')'
                )
            )
        );

//fvc::laplacian(difffield,phase1);


//dimensionedScalar valdiff("valdiff",dimdiff_,1.0e-5);
//difffield=fvc::laplacian(valdiff,phase1);




        surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];

        forAllIter(PtrDictionary<phaseModel>, phases_, iter2)
        {
            phaseModel& phase2 = iter2();
            volScalarField& alpha2 = phase2;

            if (&phase2 == &phase1) continue;

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


            surfaceScalarField phirdiff(phase1.phi() - phase1.phi());


            scalarCoeffSymmTable::const_iterator cAlpha
            (
                cAlphas_.find(interfacePair(phase1, phase2))
            );

            if (cAlpha != cAlphas_.end())
            {
                surfaceScalarField phic
                (
                    (mag(phi_) + mag(phir))/mesh_.magSf()
                );
//fvc::interpolate(epsilon1)*
                phir += (min(cAlpha()*phic, max(phic))*nHatf(phase1, phase2));
            }

//difffield=(fvc::laplacian(coefffield,phase1));


            scalarCoeffSymmTable::const_iterator dAlpha
            (
                dAlphas_.find(interfacePair(phase1, phase2))
            );

            if (dAlpha != dAlphas_.end())
            {
            dimensionedScalar valdiff("valdiff",dimdiff_,dAlpha());
                difffield=fvc::laplacian(valdiff,phase1);


            }




            word phirScheme
            (
                "div(phir," + alpha2.name() + ',' + alpha1.name() + ')'
            );

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


        }





alphalaps.set(phasei,new volScalarField(difffield));//
volScalarField& alphalap=alphalaps[phasei];



        surfaceScalarField::Boundary& alphaPhiCorrBf =
            alphaPhiCorr.boundaryFieldRef();

        // Ensure that the flux at inflow BCs is preserved
        forAll(alphaPhiCorrBf, patchi)
        {
            fvsPatchScalarField& alphaPhiCorrp = alphaPhiCorrBf[patchi];

            if (!alphaPhiCorrp.coupled())
            {
                const scalarField& phi1p = phase1.phi().boundaryField()[patchi];
                const scalarField& alpha1p = alpha1.boundaryField()[patchi];

                forAll(alphaPhiCorrp, facei)
                {
                    if (phi1p[facei] < 0)
                    {
                        alphaPhiCorrp[facei] = alpha1p[facei]*phi1p[facei];
                    }
                }
            }
        }


        MULES::limit
        (
            1.0/mesh_.time().deltaT().value(),
            geometricOneField(),
            phase1,
            phi_,
            alphaPhiCorr,
            zeroField(),
            alphalap,//alphalap
            1,//fvm::laplacian(phase1)
            0,
            true
        );

        phasei++;
    }

    MULES::limitSum(alphaPhiCorrs);
//    MULES::limitSum(alphalaps);

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

    phasei = 0;

    forAllIter(PtrDictionary<phaseModel>, phases_, iter)
    {
        phaseModel& phase1 = iter();

        surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];



        volScalarField& alphalaplac= alphalaps[phasei];



        alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase1);//

        MULES::explicitSolve
        (
            geometricOneField(),
            phase1,
            alphaPhi,
            zeroField(),
            alphalaplac//zeroField()//alphalaplac
        );

        phase1.alphaPhi() += alphaPhi;

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

        sumAlpha += phase1;

        phasei++;
    }

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

    calcAlphas();
}

All the best,

Tom

[tom_flint2012]

So I've been trying to add diffusion between selected phase pairs in multiphaseeulerfoam. I've managed to add a source term to MULES that exhibits some diffusive type behavior between selected phase pairs. However, the phase fractions of the phases is now not conserved.

Does anybody have any ideas about this? Has anybody managed to add a diffusive source term to MULES?

void Foam::multiphaseSystem::solveAlphas() //const
{
    PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
    PtrList<volScalarField> alphalaps(phases_.size());




////////////////////////////////Tried to make some fields to keep diffusion in desired phases
volScalarField difffield
    (
        IOobject
        (
            "difffield",
            mesh_.time().timeName(),
            mesh_
        ),
        mesh_,
        dimensionedScalar("difffield", dimensionSet(0, 0, -1, 0, 0), 0)
    );


    volScalarField coefffield
    (
        IOobject
        (
            "coefffield",
            mesh_.time().timeName(),
            mesh_
        ),
        mesh_,
        dimensionedScalar("coefffield", dimensionSet(0, 2, -1, 0, 0), 0)
    );


    int phasei = 0;

    forAllIter(PtrDictionary<phaseModel>, phases_, iter)
    {
        phaseModel& phase1 = iter();
        volScalarField& alpha1 = phase1;

        phase1.alphaPhi() =
            dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0);

        alphaPhiCorrs.set
        (
            phasei,
            new surfaceScalarField
            (
                "phi" + alpha1.name() + "Corr",
                fvc::flux
                (
                    phi_,
                    phase1,
                    "div(phi," + alpha1.name() + ')'
                )
            )
        );



coefffield*=0.0;
difffield*=0.0;
//difffield=(fvc::laplacian(coefffield,phase1));

//surfaceScalarField test = fvm::laplacian(alpha1);

////difffield*=0.0;
//
//coeffss.set(phasei,new volScalarField(coefffield));//
//volScalarField& coeff=coefffield[phasei];

        surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];

        forAllIter(PtrDictionary<phaseModel>, phases_, iter2)
        {
            phaseModel& phase2 = iter2();
            volScalarField& alpha2 = phase2;

            if (&phase2 == &phase1) continue;

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


            surfaceScalarField phirdiff(phase1.phi() - phase1.phi());


            scalarCoeffSymmTable::const_iterator cAlpha
            (
                cAlphas_.find(interfacePair(phase1, phase2))
            );

            if (cAlpha != cAlphas_.end())
            {
                surfaceScalarField phic
                (
                    (mag(phi_) + mag(phir))/mesh_.magSf()
                );
//fvc::interpolate(epsilon1)*
                phir += (min(cAlpha()*phic, max(phic))*nHatf(phase1, phase2));
            }

//difffield=(fvc::laplacian(coefffield,phase1));


            scalarCoeffSymmTable::const_iterator dAlpha
            (
                dAlphas_.find(interfacePair(phase1, phase2))
            );

            if (dAlpha != dAlphas_.end())
            {
            dimensionedScalar valdiff("valdiff",dimdiff_,dAlpha());
//                difffield=fvc::laplacian(valdiff,alpha1);
        coefffield=alpha1*valdiff;

            }

            word phirScheme
            (
                "div(phir," + alpha2.name() + ',' + alpha1.name() + ')'
            );

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

difffield+=(alpha2)*fvc::laplacian(coefffield,alpha1);
        }

//        difffield+=1.0e-2*alpha1*(1.0/mesh_.time().deltaT());
//rDeltaT.field();
//;

alphalaps.set(phasei,new volScalarField(difffield));//
volScalarField& alphalap=alphalaps[phasei];



        surfaceScalarField::Boundary& alphaPhiCorrBf =
            alphaPhiCorr.boundaryFieldRef();

        // Ensure that the flux at inflow BCs is preserved
        forAll(alphaPhiCorrBf, patchi)
        {
            fvsPatchScalarField& alphaPhiCorrp = alphaPhiCorrBf[patchi];

            if (!alphaPhiCorrp.coupled())
            {
                const scalarField& phi1p = phase1.phi().boundaryField()[patchi];
                const scalarField& alpha1p = alpha1.boundaryField()[patchi];

                forAll(alphaPhiCorrp, facei)
                {
                    if (phi1p[facei] < 0)
                    {
                        alphaPhiCorrp[facei] = alpha1p[facei]*phi1p[facei];
                    }
                }
            }
        }


        MULES::limit
        (
            1.0/mesh_.time().deltaT().value(),
            geometricOneField(),
            phase1,
            phi_,
            alphaPhiCorr,
            alphalap,
            alphalap,//alphalap
            1,//fvm::laplacian(phase1)
            0,
            true
        );

        phasei++;
    }

    MULES::limitSum(alphaPhiCorrs);
//    MULES::limitSum(alphalaps);

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

    phasei = 0;

    forAllIter(PtrDictionary<phaseModel>, phases_, iter)
    {
        phaseModel& phase1 = iter();

        surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];



        volScalarField& alphalaplac= alphalaps[phasei];

//alphalaplac+=;

        alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase1);//

        MULES::explicitSolve
        (
            geometricOneField(),
            phase1,
            alphaPhi,
            alphalaplac,
            alphalaplac//zeroField()//alphalaplac
        );

        phase1.alphaPhi() += alphaPhi;

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

        sumAlpha += phase1;

        phasei++;
    }

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

    calcAlphas();
}

I have the solver and test-cases for anyone who is interested.

[tom_flint2012]

To anyone reading this thread. I managed to solve this issue. I had neglected the diffusive flux terms between the offending phases.

After correcting this the volume fractions are nearly conserved and diffusion only occurs between the selected pairs of phases.

If anyone is interested in the solver, you can drop me an email at Thomas.flint@manchester.ac.uk and I’d be happy to share the solution and my test cases.

All the best,
Tom

[wavefunction]

If there were only 3 phases being used, you could have tried using interMixingFoam as well as the solver allows diffusion between two phases.

[wavefunction]

I am having trouble adding a source term in MULES. I want to modify interMixingFoam to accomodate phase change. So I want to add the source terms -mDot/rho1 to alpha1Eqn and mDot/rho2 to alpha2Eqn.

mDot = D*rhogp*mag(fvc::grad(Y))*mag(fvc::grad(alpha1))/(scalar(1) - Y)

where D - diffusion coefficient, Y - mass fraction, rhogp - gas phase density.

Without adding the source terms to MULES, my case using my custom solver is diverging immensely. On running gdb, I found that the problem lies at line 442 of MULESTemplates.C.

Please help me.

Program received signal SIGFPE, Arithmetic exception.
0x00000000004900b1 in Foam::MULES::limiter<double, Foam::geometricOneField, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>, Foam::zeroField> (
allLambda=..., rDeltaT=@0x7fffffff4cb0: 400000000, rho=..., psi=...,
phiBD=..., phiCorr=..., Sp=..., Su=..., psiMax=1, psiMin=0)
at /home/saicharanb56/OpenFOAM/OpenFOAM-5.x/src/finiteVolume/lnInclude/MULESTemplates.C:442
442 max(min

[tom_flint2012]

Late response but if you are still battling this you need to think about your Su and Sp fields in mules and then correct the relative fluxes.