CFDEM Particle Forces and Their Implementation in the Code.

Particle forces in CFDEM are given by

{m_p}{{d{V_p}} \over {dt}} = {m_p}g + {f_{\nabla \cdot \tau }} + {f_{\nabla p}} + \sum\limits_{{N_p}} {{f_{p - p}} + } \sum\limits_{{N_p}} {{f_{p - w}}} + {f_{drag}}

Particle Volume / Void Fraction

It is first prudent to discuss how the volume of the particle is calculated in CFDEM because this volume is used in almost all particle forces.

{V_s} is the volume of partial or full particle (depending on void fraction model) volume. Looking at {V_s} in the code we see that it is calculated by 

Vs = particleCloud_.particleVolume(index)

Where {V_s} is the particle partial or full volume. Now this will depend on the void fraction model you used and so looking at {V_s} in the code we see that it is calculated by 

Vs = particleCloud_.particleVolume(index)

The particleVolume function is calculated by 

inline double cfdemCloud::particleVolume(int index) const
{
    return particleV_[index][0];
}

In #include "cfdemCloudI.H". Where cfdemCloudI.H is included by cfdemCloud.H at the end of it’s code.

If we go through the code we have starting at gradPForce.H -> forceModel.H-> cfdemCloud.H ->"cfdemCloudI.H" to end up at particleVolume().

Now that we have particleVolume() we know it returns particleV_. We find particleV_ in dividedVoidFraction.C and setWeightedSource.H.

This is where things mildy interesting. We see that in dividedVoidFraction.C there is a loop starting at 

for (int index=0; index < particleCloud_.numberOfParticles(); index++)

All particleWeights, particleVolumes, and particleV are set equal to 0 in this loop initially and so every time this loops it is reset and the new particle partial/full volumes are calculated. 

cellID = particleCloud_.cellIDs()[index][0];

Particle volume inside of the loop is calculated by volume = Vp(index,radius,scaleVol); with Vp() located in dividedVoidFraction.H as 

       return constant::mathematical::fourPiByThree*radius*radius*radius*scaleVol;

cellID associated with that particle index is found by 

cellID = particleCloud_.cellIDs()[index][0];

voidFractionNext_ is then calculated by using cellID associated with that particle index using 

                scalar newAlpha = voidfractionNext_[cellID]- volume*centreWeight/cellVol;
                if (newAlpha > alphaMin_)
                {
                    voidfractionNext_[cellID] = newAlpha;
                }
                else
                {
                    voidfractionNext_[cellID] = alphaMin_;
                    tooMuch_ += (alphaMin_-newAlpha) * cellVol;
                }

Finally, particleVolumes and particleV is stored by

                particleVolumes[index][0] += volume*centreWeight;
                particleV[index][0] += volume*centreWeight;

Where particleV is the return of the function particleVolume().

Pressure Force, {f_{\nabla p}} = - {V_s}\nabla p

The pressure force is given by {f_{\nabla p}} = - {V_s}\nabla p in gradPForce.C and gradPForce.H

This force is relatively straight forward. In gradPForce.C the pressure is set to

p_(sm.mesh().lookupObject<volScalarField> (pFieldName_))

Then gradient is calculated of the pressure in the respective cell and multiplied by the particle volume, Vs, to get the force. It should be noted that there is a pressure interpolation function in CFDEM but I won’t go into this in the post.

volVectorField gradPField = fvc::grad(p_);
gradP = gradPField[cellI]
force = -Vs*gradP*(1.0+addedMassCoeff_);

Where addedMassCoeff is an artificial mass added by user in pressure properties of couplingProperties file.

Viscous Force

The viscous force is given by

\tau = {\rho _f}\nu \left( {\nabla u + {{\left( {\nabla u} \right)}^T}} \right)

Starting at viscForce.C we see that the viscous force is calculated by

force = -Vs*divTau*(1.0+addedMassCoeff_);

Calculation of divTau is located in the same file with 

const volVectorField& divTauField = forceSubM(0).divTauField(U_);
divTau = divTauField[cellI]

We then find the divTauField() function in forceSubModel.C

const volVectorField& forceSubModel::divTauField(const volVectorField& U) const
{
    // calc div(Tau)
    #ifdef compre
        const volScalarField& mu_ = muField();
        divTau_ = -fvc::laplacian(mu_, U) - fvc::div(mu_*dev(fvc::grad(U)().T()));
        return divTau_;
    #else
        const volScalarField& nu_ = nuField();
        const volScalarField& rho_ = rhoField();
        divTau_ = -fvc::laplacian(nu_*rho_, U)- fvc::div(nu_*rho_*dev(fvc::grad(U)().T()));
        return divTau_;
    #endif
}

where 

divTau_ = -fvc::laplacian(nu_*rho_, U)- fvc::div(nu_*rho_*dev(fvc::grad(U)().T()));

We can see that we take the Laplacian of the velocity first multiplied by the density and viscosity. Then we take the divergence of the deviatoric particle of the transform of the gradient of velocity. To put it in math terms…

{f_{\nabla \cdot \tau }} = - \nabla \cdot \left( {\nabla u + {{\left( {\nabla u} \right)}^T}} \right){\rho _f}\nu = - {\rho _f}\nu \left( {\Delta u} \right) - {\rho _f}\nu \left( {\nabla \cdot {{\left( {\nabla u} \right)}^T}} \right)

Archimedes lift force

The Archimedes lift force is another easy force that is implemented given by

{f_l} = {\rho _f}g{V_p}

The force is calculated by Archimedes’ principle that a solid body in a fluid will experience opposing force to gravity because of the difference in the density.

The force is implemented in Archimedes.C by

force = -g_.value()*forceSubM(0).rhoField()[cellI]*particleCloud_.particleVolume(index);

Drag Force and Contact Forces

This topic is quite extensive and so I’ll talk about this either by editing this post or adding a new post later on.

Leave a comment