CFDEM Coupling

We already went over the DEM forces and how they are implemented in the CFDEM code by default. I will now go over the coupling source term and how it’s implemented in the code by default. It is always good to see if their are some disparities between what they say their formulas use and what is actually used in the code. I’ve come to realize that this is critical to do in OpenSource code and is actually encouraged because, as being a human, we tend to make mistakes unfortunately. Looking into code can’t be done in most commercial software and you just have to take their word for it but when we have the source code we have no excuses!

According to the CFDEM website: https://www.cfdem.com/media/CFDEM/docu/CFDEMcoupling_Manual.html. If we quote them exactly, “‘modelType’ refers to the formulation of the equations to be solved. Choose “A”, “B” or “Bfull”, according to Zhou et al. (2010): “Discrete particle simulation of particle-fluid flow: model formulations and their applicability”, JFM. “A” requires the use of the force models gradPForce and viscForce, whereas “B” requires the force model “Archimedes”. “Bfull” refers to model type I, “A” refers to model type II and “B” refers to type III in the nomenclature used by Zhou et al.

So I take it from this sentence that the coupling of the code and that the equations stated by Zhou et al. are in the code. Taking a snapshot of the equations at the paper Zhou et al. [1] we first see mass conservation given by

This is straight forward but now we have the three sets of momentum equations given by

Set 1 taken directly from Zhou et al. 2010 (Model “Bfull” in CFDEM)

Set 2 taken directly from Zhou et al. 2010 (Model “A” in CFDEM)

Set 3 taken directly from Zhou et al. 2010 (Model “Bfull” in CFDEM)

Comparing Set 1 and Set 2 we can see right away that the momentum source term for Set 1 includes all of the DEM forces; that is the viscous, pressure, drag, and any additional forces such as a Saffman lift force where Set 2 does not. Instead, Set 2 has that the pressure gradient and viscous forces separated from the fluid-particle interaction source term and the pressure gradient and viscous forces are multiplied by volume fraction. Essentially, or shall I say mathematically, these two sets are equal but it is merely the order in which the contribution of viscous and pressure force contribution in the momentum equation is calculated. Is the pressure gradient and velocity force contribution calculated with the U equation or before the U equation with the source momentum term? Does this mean that one is explicit and one is implicit? Well lets look at the code and find out what they have done there.

I’m going to be going over the code as implimented in https://github.com/ParticulateFlow/CFDEMcoupling-PFM instead of CFDcoupling Public.

Starting in cfdemSolverPiso solver we go to the momentum equation to determine what is being solved and in what order. In Ugn.H we see that the momentum equation is solved by

fvVectorMatrix UEqn
(
    fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
  + fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
  + particleCloud.divVoidfractionTau(U, voidfraction)
  ==
    fvOptions(U)
  - fvm::Sp(Ksl/rho,U)
);

In through term by term..

${\rm{fvm::ddt(voidfraction,U) }} \to {{\partial \left( {{\varepsilon _f}u} \right)} \over {\partial t}}$

${\rm{fvm::div(phi,U)}} \to \nabla \cdot \left( {{\varepsilon _f}uu} \right)$

where phi is calculated by voidfraction*U in createFields.H

Now looking at ${\rm{particleCloud}}{\rm{.divVoidfractionTau(U, voidfraction)}}$ we see that this function is in cfdemCloud.C with the return of

    return
    (
      - fvm::laplacian(voidfractionNuEff(voidfraction), U)
      - fvc::div(voidfractionNuEff(voidfraction)*dev2(fvc::grad(U().T()))
    );

This is the expanded out version of the divergence of the fluid stress tensor.

${\varepsilon _f}\nabla \cdot \tau = - {\varepsilon _f}\nabla \cdot \left( {\nabla u + {{\left( {\nabla u} \right)}^T}} \right) = - {\varepsilon _f}\left( {\Delta u} \right) - {\varepsilon _f}\left( {\nabla \cdot {{\left( {\nabla u} \right)}^T}} \right) \to$

After the specification we see an if statement to switch between model “Bfull” and “A” or “B” that is with

if (piso.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
{
    solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
    fvOptions.correct(U);
}
else if (piso.momentumPredictor())
{
    solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
    fvOptions.correct(U);
}

It should first be noted that there is the use of the OpenFOAM function Sp() with terms being added to the LHS and RHS of the momentum equation implicitly to contribute to the diagonal of the matrix and help with convergence. This is done because we have explicit source terms added from the DEM side that could make our matrix ill-posed.

Tying it all together in math form the full equation is given by

${{\partial \left( {{\varepsilon _f}{u_l}} \right)} \over {\partial t}} + \nabla \cdot \left( {{\varepsilon _f}{u_l}{u_l}} \right) = - switch \cdot \nabla p - Ksl\left( {{u_l} - {u_s}} \right) + \nabla \cdot \left( {switch \cdot \tau } \right) + f$

where $switch$ is is equal to void fraction, ${{\varepsilon _f}}$ , for model A (set II) and $1$ for model Bfull (set II).

Mathematically the model A and model Bfull are equivalent but it does change the order in which particle to fluid momentum coupling occurs. Model A can be seen as implicit in that it goes against Newton’s third law by taking in the volumetric fraction of solid momentum while calculating the fluid velocity where Model Bfull uses the fluid velocity from the previous time step to calculate the coupling source term. This being said Zhou et. al 2010 does state that there is very little difference between the two models.

Now we will go into the calculation of the momentum exchange term $Ksl$ in the code. In CFDEM the calculation of $Ksl$ is given by

$Ksl = {{{\varepsilon _f}\left| {\sum\limits_i {{F_d}} } \right|} \over {{V_c}\left| {{u_l} - {u_p}} \right|}}$

where ${\sum\limits_i {{F_d}} }$ is the summation of all the particle to fluid interaction forces. To clarify we have for Model A (set II)

$\sum\limits_i {{F_d}} = {f_d} + {f^{''}}$

where ${f_d}$ is the particle drag force and ${f^{''}}$ is the summation of all other forces besides viscous and pressure forces. For example a Saffman lift force would be used here.

For Model Bfull (set I) we have

References

$\sum\limits_i {{F_d}} = {f_d} + {f_{\nabla \tau }} + {f_{\nabla p}} + {f^{''}}$

where ${f_{\nabla \tau }}$ is the viscous force due to intermolecular forces between the particles and fluid, and ${f_{\nabla p}}$ is the pressure gradient force.

Looking at where $Ksl$ in the code we see that first the particles are “evolved” or moved on the DEM side in cfdemSolverPiso.C by

        particleCloud.clockM().start(2,"Coupling");
        bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);

Once particles are evolved then the field is first smoothed and $Ksl$ is calculated by

Ksl = particleCloud.momCoupleM(0).impMomSource();

We see that the impMomSource() function is called. This function is called in implicitCoupling.C with

tmp<volScalarField> implicitCouple::impMomSource()

Ksl in then calculated by

    return tmp<volScalarField>
    (
        new volScalarField("Ksl_implicitCouple", (1. - tsf)*KslPrev_ + tsf * KslNext_)
    );

Where KslNext is calculated previously in the same file by

KslNext_[cellI] = mag(particleCloud_.forceM(0).impParticleForces()[cellI])/Ur/particleCloud_.mesh().V()[cellI];

KslPrev is set when resetMomSourceField() is called in cfdemCloud.C with the function given by

    return tmp<volScalarField>
    (
        new volScalarField("Ksl_implicitCouple", (1. - tsf) * KslPrev_ + tsf * KslNext_)
    );

tsf is calculated by

const scalar tsf = particleCloud_.dataExchangeM().timeStepFraction();

with timeStepFraction() function called in dataExchangeModel.C for returning the fraction between previous coupling time step and actual time step. This ensures that Ksl is always an interpolated value from the previous coupling time step and the current time step.

Now back to the Ksl term. We can see that KslNext_[cellI] is calculated with impParticleForces(). This is set in cfdemCloud.C by adding all the impForces_ inside of each cell by

    averagingM().setVectorSum
    (
        forceM(0).impParticleForces(),
        impForces_,
        particleWeights_,
        NULL //mask
    );

Ok, so now where/how are impForces calculated? We see in forceSubModel.C that it is calculated in the member fucntion partToArray() by

    if(!switches_[SW_TREAT_FORCE_DEM])// !treatDEM
    {
        if(switches_[SW_TREAT_FORCE_EXPLICIT]) // treatExplicit
        {
            for(int j=0;j<3;j++)
                particleCloud_.expForces()[index][j] += dragTot[j];
        }
        else   //implicit treatment, taking explicit force contribution into account
        {
            for(int j=0;j<3;j++)
            {
                particleCloud_.impForces()[index][j] += dragTot[j] - dragEx[j]; //only consider implicit part!
                particleCloud_.expForces()[index][j] += dragEx[j];
            }
        }
    }

Essentially this uses a switch to determine if coupling force for viscosity and pressure are used. If treatForceExplicit is specified in couplingProperties expForces will only be used with adding all of the coupling forces together. Otherwise, the the explicit drag calculation will be subtracted out where dragEx is the third term in the function partToArray. For example if we look into how drag is calculated in the GidaSpowDrag.C we can see that both drag and an explicit drag is calculated.

forceSubM(0).partToArray(index,drag,dragExplicit,Ufluid,dragCoefficient);

So now, when is this switch turned on? It appears this doesn’t matter anyway because impForces() is always used by default in CFDEM. Therefore we need to look at how the explicit forces using dragEx

Therefore looking back to GidaspowDrag.C we can see an explicit correction that sets the drag to explicit given by

forceSubM(0).explicitCorr(drag,dragExplicit,dragCoefficient,Ufluid,U_[cellI],Us,UsField_[cellI],forceSubM(0).verbose());

Now we have the explicitCorr function because this is where dragExplicit is calculated and ultimately whether we subtract that portion from impForces or not. The explicitCorr function is located in forceSubModel.C with it given by

void forceSubModel::explicitCorr
(
    vector& dragImplicit,
    vector& dragExplicit,
    scalar& dragCoefficient,
    vector& Ufluid,
    const vector& Ucell,
    vector& Us,
    const vector& UsCell,
    bool verbose,
    label index
) const
{
    dragExplicit=vector::zero;
}

So essentially this function merely sets dragExplicit to zero and therefore it appears that by default in CFDEM we cannot use any specification of explicit forces or rather use Bfull model for that matter. Now, I keep talking about drag force which might be confusing because this is actually the total coupling forces used in the Ksl term because the partToArray function is used at the end of every force call. Lets look at the pressure and viscous forces to illustrate this. In viscForce.C and gradPForce.C we see at the end of the code that the force is added with

forceSubM(0).partToArray(index,force,vector::zero);

In conclusion we have that as each force is calculated an “if” statement is used inside of each force file (gradPForce.C, GidaspowDrag.C,viscForce.C...) to set the SW_TREAT_FORCE_DEM switch to true or false. An if statement is then used with this switch in forceSubModel.C to determine if the force is added to the impForces() which ultimately….. -> impMomSource() ->Ksl.

We see that code is indeed following the equations as outlined in the Zhou et al. paper. This is good news and I can precede. It’s always good know that the equations think are being solved are actually being solved. Especially when there is a free and full (paid) version to software. You don’t know if some features are turned off or wrong and that the full version might be correct.

[1] Zhou, Z. Y., Kuang, S. B., Chu, K. W., & Yu, A. B. (2010). Discrete particle simulation of particle–fluid flow: Model formulations and their applicability. Journal of Fluid Mechanics, 661, 482–510. https://doi.org/10.1017/S002211201000306X

Share this:

Related

Leave a comment Cancel reply