[keepfit]

Dear Foamers,

Having been working on solid particles - free surface flow interaction for a while, so far I got some good results in terms of coupled particle-fluid interaction. Nevertheless, an unnoticeable issue was found on the volume conservation of fluid phase. The modified solver is based on interFoam where an open-source DEM (discrete element method) code is called as a shared library to handle dynamics of solid particles. The interaction is realized via a momentum exchange term (Sp_f) to the R.H.S of N-S equation, and empirical drag models to calculate fluid forces acting on solid particles.

The VOF and continuity equations used in the solver are given by:

hence the Source term which stands for the displaced volume by solid particles reads as:


here "alpha_f" is the volume fraction of fluid phase in a cell (= 1- V_particles/V_cell), and "U_f" the fluid velocity. The continuity equation is already implemented into the pressure equation instead of the original one.

The alplaEqn.H is implemented as:

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

    surfaceScalarField phic(mag(phi/mesh.magSf()));
    phic = min(interface.cAlpha()*phic, max(phic));
    surfaceScalarField phir(phic*interface.nHatf());

    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
    {   
        surfaceScalarField phiAlpha  
        (
            fvc::flux
            (
                phi,  // = U*alpha_f
                alpha1,
                alphaScheme
            )
          + fvc::flux
            (
                -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
                alpha1,
                alpharScheme
            )
        );

       //MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);  

       // source term
       volScalarField divU(fvc::div(U*alpha_f));
       volScalarField SpVoF = alpha1*divU;

       MULES::explicitSolve(geometricOneField(), alpha1, phi, phiAlpha, SpVoF, zeroField(), 1, 0);

       rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
    }  // end for

    Info<< "Phase-1 volume fraction = "
        << alpha1.weightedAverage(mesh.Vsc()).value()
        << "  Min(alpha1) = " << min(alpha1).value()
        << "  Max(alpha1) = " << max(alpha1).value()
        << endl;
}

The solver runs well in general, the behavior of solid particles under the influence of fluid phase in the simulation is expected. The following test is to verify the modified VOF method on the volume conservation: a number of solid particles was poured into a box with half-filled water (0.1m high), after the particles rest on the bottom, the water level will rise.


However. what i got from the coupled particle-fluid simulation is odd: the water level remains at 0.1 m which is not true after several hundred of particles poured into the water column.


Where is the problem in the VOF implementation, especially the source term and the MULES solver. I read the MULES.c and MULESTemplates.C but still could not fully understand.

Can we directly implement the alpha equation just like in UEqn.H without using MULES (e.g. bubbleFoam)? I would appreciate any helpful advice

Cheers,
David

[keepfit]

Update

The direct implementation of the modified VOF method seems working, but still suffer some mass loss.

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

    surfaceScalarField phic(mag(phi/mesh.magSf()));
    phic = min(interface.cAlpha()*phic, max(phic));
    surfaceScalarField phir(phic*interface.nHatf());
    
    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
    {

        volScalarField::DimensionedInternalField SpVof
        (
            IOobject
            (
                "SpVof",
                runTime.timeName(),
                mesh
            ),
            // Divergence term is handled explicitly to be
            // consistent with the explicit transport solution
            fvc::div(alpha_f*U) 
        );

        fvScalarMatrix alpha1Eqn
        (
             fvm::ddt(alpha1)
           + fvm::div(phi, alpha1, alphaScheme)
           + fvm::div(-fvc::flux(-phir, alpha_f - alpha1, alpharScheme), alpha1, alpharScheme)  // alpha1 + alpha2 = alpha_f;  alpha_f +alpha_solid = 1
           == fvm::Sp(SpVof, alpha1)
        );
        alpha1Eqn.relax();
        alpha1Eqn.solve();

     rhoPhi = alpha1Eqn.flux()*(rho1 - rho2) + phi*rho2;
    }  // end for

    Info<< "Phase-1 volume fraction = "
        << alpha1.weightedAverage(mesh.Vsc()).value()
        << "  Min(alpha1) = " << min(alpha1).value()
        << "  Max(alpha1) = " << max(alpha1).value()
        << endl;
}

any advice?

thanks!

[ramakant]

Dear David Long

this model is also used for the lake modelin ?.For example in the lake water and solid particle are coming together and get settled at the bottom of the lake/river.

any advice?
thanks!

[keepfit]

Originally Posted by ramakant

Dear David Long

this model is also used for the lake modelin ?.For example in the lake water and solid particle are coming together and get settled at the bottom of the lake/river.

any advice?
thanks!

yes, the solver can be used to study river/channel sediment transport as well.

[Dan M]

Hello,

I know this is an old thread but I'm dealing with exactly the same problem with water/air/particles solver (i.e. no volume displacement due to addition of particles). I'm able to maintain the volume fraction of the liquid phase and keep alpha1 between 0 and 1, but the water level doesn't rise accordingly. I've tried MULES and direct implementation.

Have you managed back then to implement the alpha equation correctly? Could you share your results?

Your help would be much appreciated. Thanks!

[gridley2]

Hi, did you try adding your alpha source term in alphaSuSp.H? That will make it get used by MULES correctly. There are a few equations where consistent source terms need to be added in for MULES, which is why they seem to have added that alphaSuSp file.