[Friendly]

Hi,

I try to understand how OpenFOAM calculates the drag forces and coefficients. In the sourceCode of forceCoeffs.C I have found the following exression:

bool Foam::functionObjects::forceCoeffs::write()
{
    forces::calcForcesMoment();

    if (Pstream::master())
    {
        writeFiles::write();

        scalar pDyn = 0.5*rhoRef_*magUInf_*magUInf_;

        Field<vector> totForce(force_[0] + force_[1] + force_[2]);
        Field<vector> totMoment(moment_[0] + moment_[1] + moment_[2]);

        List<Field<scalar>> coeffs(3);
        coeffs[0].setSize(nBin_);
        coeffs[1].setSize(nBin_);
        coeffs[2].setSize(nBin_);

        // lift, drag and moment
        coeffs[0] = (totForce & liftDir_)/(Aref_*pDyn);
        coeffs[1] = (totForce & dragDir_)/(Aref_*pDyn);
        coeffs[2] = (totMoment & pitchAxis_)/(Aref_*lRef_*pDyn);

        scalar Cl = sum(coeffs[0]);
        scalar Cd = sum(coeffs[1]);
        scalar Cm = sum(coeffs[2]);

Nothing special here, everything is more or less as expected. I just need to find out how the force [0],[1] and [2] are calculated. I have found out that 0 is the normal, 1 the tangential and 2 the porous force. But what is a porous force? I have never read such thing before.

Looking in forces.C I have found:

if (directForceDensity_)
    {
        const volVectorField& fD = obr_.lookupObject<volVectorField>(fDName_);

        const fvMesh& mesh = fD.mesh();

        const surfaceVectorField::Boundary& Sfb =
            mesh.Sf().boundaryField();

        forAllConstIter(labelHashSet, patchSet_, iter)
        {
            label patchi = iter.key();

            vectorField Md
            (
                mesh.C().boundaryField()[patchi] - coordSys_.origin()
            );

            scalarField sA(mag(Sfb[patchi]));

            // Normal force = surfaceUnitNormal*(surfaceNormal & forceDensity)
            vectorField fN
            (
                Sfb[patchi]/sA
               *(
                    Sfb[patchi] & fD.boundaryField()[patchi]
                )
            );

            // Tangential force (total force minus normal fN)
            vectorField fT(sA*fD.boundaryField()[patchi] - fN);

            //- Porous force
            vectorField fP(Md.size(), Zero);

            applyBins(Md, fN, fT, fP, mesh.C().boundaryField()[patchi]);
        }
    }

So forces are calculated via a forceDensity, but how does the calculation of the forceDensity work?

Does anyone know how OpenFOAM calculates the forces and coefficients? I would be very thankfull for any help.


Greetings,
Friendly