CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Add diffusion between selected phase pairs in multiphaseeulerfoam

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By tom_flint2012

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   March 12, 2018, 10:16
Default Add diffusion between selected phase pairs in multiphaseeulerfoam
  #1
Member
 
Thomas Flint
Join Date: Jan 2016
Posts: 56
Rep Power: 9
tom_flint2012 is on a distinguished road
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!


Code:
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 is offline   Reply With Quote

Old   March 14, 2018, 09:33
Default Update:A step in the right direction?
  #2
Member
 
Thomas Flint
Join Date: Jan 2016
Posts: 56
Rep Power: 9
tom_flint2012 is on a distinguished road
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?

Code:
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.
Attached Images
File Type: jpg diffusion_multiphaseEulerFoam.jpg (25.4 KB, 43 views)
Stefanie.S.W. likes this.
tom_flint2012 is offline   Reply With Quote

Old   March 15, 2018, 16:49
Default Resolved
  #3
Member
 
Thomas Flint
Join Date: Jan 2016
Posts: 56
Rep Power: 9
tom_flint2012 is on a distinguished road
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 Id be happy to share the solution and my test cases.

All the best,
Tom
tom_flint2012 is offline   Reply With Quote

Old   March 20, 2018, 08:41
Default interMixingFoam
  #4
New Member
 
Saicharan
Join Date: Jan 2018
Location: Bangalore, India
Posts: 29
Rep Power: 7
wavefunction is on a distinguished road
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 is offline   Reply With Quote

Old   March 20, 2018, 08:45
Default Adding Source term in MULES
  #5
New Member
 
Saicharan
Join Date: Jan 2018
Location: Bangalore, India
Posts: 29
Rep Power: 7
wavefunction is on a distinguished road
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.

Quote:
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
wavefunction is offline   Reply With Quote

Old   February 6, 2019, 07:11
Default MULES source terms
  #6
Member
 
Thomas Flint
Join Date: Jan 2016
Posts: 56
Rep Power: 9
tom_flint2012 is on a distinguished road
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.
tom_flint2012 is offline   Reply With Quote

Reply

Tags
diffusion, interface, mules, multiphase, multiphaseeulerfoam

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
InterFoam: Add an equation to only one phase voingiappone OpenFOAM Programming & Development 41 June 7, 2022 09:04
[PyFoam] and paraview eelcovv OpenFOAM Community Contributions 28 May 30, 2016 09:23
diffusion Term for add. Variable Zaktatir CFX 9 July 16, 2011 09:51
Source Term used in Eulerian Model(Two phase) Padian FLUENT 1 May 19, 2008 03:47
add user scalar in one phase zhu CFX 0 April 27, 2002 03:45


All times are GMT -4. The time now is 15:31.