need some help for adding polydisperse drag force

[Nimajhi]

I want to add the Beetstra drag force for a polydisperse system. Implementing the drag force is easy, but the problem is with the "Mean Sauter Diameter" that needs to be in the Reynolds calculation. For the Mean Sauter Diameter, we need to access the solid volume fraction of each particle and the total of them, and the diameter of each particle.
Please give some guidance in this regard

[PhasicFlow]

It requires you to write have access to the diameter of particles only and you can loop over particles in each cell to obtain the mean diameter. You have a map from particle index to cell index, while for this calculation you need a map from cell index to all particles inside that cell. You have two options:

1. Construct this map and calculate the mean diameter for each cell and then rest of calculations.

2. Construct a vector with size of the cells, and loop over all cells and to see which particle is in which cell and how it contributes to the mean diameter. I guess you need to perform loop over particles twice here.

Give it a try and comeback to me if you needed any further help.

[Nimajhi]

I doubt this code is completely correct, please take a look at it

#include "sauterMeanDiameter.hpp"
#include "couplingMesh.hpp"
#include "streams.hpp"

namespace pFlow::coupling
{

sauterMeanDiameter::sauterMeanDiameter(const couplingMesh& cMesh)
:
    cMesh_(cMesh),
    cellToParticles_(cMesh.nCells()),
    d32_(cMesh.nCells(), 0.0)
{}

void sauterMeanDiameter::calculate
(
    const labelList& parCells,
    const Plus::realProcCMField& diameter
)
{
    //Reset

    #pragma omp parallel for
    for (label c = 0; c < cMesh_.nCells(); ++c)
    {
        cellToParticles_[c].resize(0); 
        d32_[c] = 0.0;
    }

    //cell to particle map

    for (label p = 0; p < parCells.size(); ++p)
    {
        const label c = parCells[p];
        if (c >= 0 && c < cMesh_.nCells())
        {
            cellToParticles_[c].append(p);
        }
    }

    //Calculate Sauter mean diameter for each cell

    #pragma omp parallel for
    for (label c = 0; c < cMesh_.nCells(); ++c)
    {
        const auto& particles = cellToParticles_[c];
        const label nParticles = particles.size();

        if (nParticles == 0)
        {
            d32_[c] = 0.0;
            continue;
        }

        scalar numSum = 0.0; // Sum of d^3
        scalar denSum = 0.0; // Sum of d^2

        for (label i = 0; i < nParticles; ++i)
        {
            const label p = particles[i];
            const scalar d = diameter[p];
            numSum += d * d * d;  //
            denSum += d * d;      // 
        }

        d32_[c] = (denSum > SMALL) ? (numSum / denSum) : 0.0;
    }
}

const List<DynamicList<label>>& sauterMeanDiameter::cellToParticles() const
{
    return cellToParticles_;
}

const scalarField& sauterMeanDiameter::d32() const
{
    return d32_;
}

} // pFlow::coupling

[PhasicFlow]

I think this code works just well. But you must also check the results by some tests. Create an ordered array of particles with similar size and also with two different sizes (consider every odd particle small and every even particle large) and then calculate the mean diameter.

If it was successful, then keep it.
Later you may want to check the second method i explained is more efficient or not (especially if this step is very time consuming).

Go ahead and write the code for the rest of calculations.