[sven]

Hi,

I am simulating a turbulent flow around a square. For This flow, I want to calculate the lift and the drag coefficients CD and CL. By scanning through posts in this forum, I found that this is possible by adding the following lines to the controlDict file:

functions
(
forces
{
type forces;
functionObjectLibs ("libforces.so"); //Lib to load
patches (cylinder); // change to your patch name
rhoInf 1.204; //Reference density for fluid
CofR (0 0 0); //Origin for moment calculations
outputControl timeStep;
outputInterval 100;
}

forceCoeffs
{
type forceCoeffs;
functionObjectLibs ("libforces.so");
patches (cylinder); //change to your patch name
rhoInf 1.204;
CofR (0 0 0);
liftDir (0 1 0);
dragDir (1 0 0);
pitchAxis (0 0 0);
magUInf 0.1;
lRef 1;
Aref 1;
outputControl timeStep;
outputInterval 100;
}
);

Unfortunately this does not yield the expected values for CL and CD. For fixing this problem and obtaining the correct CL/CD I first need some information on the code posted above:

1. rhoInf: What is the reference density? The density in my free stream?
2. What is magUInf? How can i determine/calculate this value?
3. The program seems to use "libforces.so". I could find this file, but where can I find the corresponding source code?
4. Where are the equations for calculating CD and CL defined?

I appreciate your help. Thank you very much

[Simon Lapointe]

Hi,

1. rhoInf is the freestream density.
2. magUinf is the magnitude of the freestream velocity.
3. The source code for the forces and forceCoeffs are located, respectively, in OpenFOAM-1.6/src/postProcessing/forces/ and OpenFOAM-1.6/src/postProcessing/forceCoeffs/
4. The equations of CL and CD are defined in OpenFOAM-1.6/src/postProcessing/forceCoeffs/forceCoeffs.C

Hope this helps

[sven]

Thanks for your answer Simon. I still have a problem with magUInf, the magnitude of the free stream velocity. What exactly is the magnitude of my free stream velocity? My free stream velocity is 0.322, so the magnitude is 0.1? I tried it with that value (0.1), but it does not yield the correct results! Thanks again for your help!

[bigred]

If your free stream velocity is 0.322, then I'm guessing the magnitude is 0.322. If that is the only velocity component ie, value (0.322 0 0). Magnitude is calculated using typical vector mathematics v²=vx²+vy²+vz². The square of the velocity magnitude is the sum of the squares of your velocity components.

[sven]

Ok..now I got it! Thank you all very much for your help. Greetings from California

Sven

[sven]

Well, I looked up the definition for CD and CL in the source Code. First of all, the right location for the files is:
/home/user/OpenFOAM/OpenFOAM-1.6/src/postProcessing/functionObjects/forces

There are two folders (forces and forceCoeffs) where you can find the appropriate .C and .H files.

I first had a look at the forceCoeffs.C file, here is part of that:

void Foam::forceCoeffs::write()
{
    if (active_)
    {
        // Create the forces file if not already created
        makeFile();

        forcesMoments fm = forces::calcForcesMoment();

        scalar pDyn = 0.5*rhoRef_*magUInf_*magUInf_;

        vector totForce = fm.first().first() + fm.first().second();
        vector totMoment = fm.second().first() + fm.second().second();

        scalar liftForce = totForce & liftDir_;
        scalar dragForce = totForce & dragDir_;
        scalar pitchMoment = totMoment & pitchAxis_;

        scalar Cl = liftForce/(Aref_*pDyn);
        scalar Cd = dragForce/(Aref_*pDyn);
        scalar Cm = pitchMoment/(Aref_*lRef_*pDyn);

        if (Pstream::master())
        {
            forcesFilePtr_()
                << obr_.time().value() << tab
                << Cd << tab << Cl << tab << Cm << endl;

            if (log_)
            {
                Info<< "forceCoeffs output:" << nl
                    << "    Cd = " << Cd << nl
                    << "    Cl = " << Cl << nl
                    << "    Cm = " << Cm << nl
                    << endl;
            }
        }
    }
}

OpenFOAM first calculates a total force, by multiplying it with the directions of each component, this force is divided into its components (lift/drag). These components are then divided by the referenceDensity multiplied with the refArea which leads to CL and CD.
Unfortunately I do not understand, how OpenFOAM calculates the total force mentioned above. I tried to look up the definition of first() and second(), which are defined in Tuple2.H (OpenFoam/InInclude/Tuple2.H). Here is the definition of these functions:

//- Return first
        inline const Type1& first() const
        {
            return f_;
        }

        //- Return first
        inline Type1& first()
        {
            return f_;
        }

        //- Return second
        inline const Type2& second() const
        {
            return s_;
        }

        //- Return second
        inline Type2& second()
        {
            return s_;
        }

I don't understand the definitions of these functions. Can anyone help me with that?

[sven]

In the forces.C file, I also found this:

Foam::forces::forcesMoments Foam::forces::calcForcesMoment() const
{
    forcesMoments fm
    (
        pressureViscous(vector::zero, vector::zero),
        pressureViscous(vector::zero, vector::zero)
    );

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

        const fvMesh& mesh = fD.mesh();

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

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

            vectorField Md = mesh.C().boundaryField()[patchi] - CofR_;

            scalarField sA = mag(Sfb[patchi]);

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

            fm.first().first() += sum(fN);
            fm.second().first() += sum(Md ^ fN);

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

            fm.first().second() += sum(fT);
            fm.second().second() += sum(Md ^ fT);
        }
    }
    else
    {
        const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
        const volScalarField& p = obr_.lookupObject<volScalarField>(pName_);

        const fvMesh& mesh = U.mesh();

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

        tmp<volSymmTensorField> tdevRhoReff = devRhoReff();
        const volSymmTensorField::GeometricBoundaryField& devRhoReffb
            = tdevRhoReff().boundaryField();

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

            vectorField Md = mesh.C().boundaryField()[patchi] - CofR_;

            vectorField pf = Sfb[patchi]*p.boundaryField()[patchi];

            fm.first().first() += rho(p)*sum(pf);
            fm.second().first() += rho(p)*sum(Md ^ pf);

            vectorField vf = Sfb[patchi] & devRhoReffb[patchi];

            fm.first().second() += sum(vf);
            fm.second().second() += sum(Md ^ vf);
        }
    }

    reduce(fm, sumOp());

    return fm;
}

Can anyone by this understand, how OpenFOAM calculates these forces?

Thanks a lot for all your help!

[desert_1250]

Originally Posted by sven

In the forces.C file, I also found this:

Foam::forces::forcesMoments Foam::forces::calcForcesMoment() const
{
    forcesMoments fm
    (
        pressureViscous(vector::zero, vector::zero),
        pressureViscous(vector::zero, vector::zero)
    );

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

        const fvMesh& mesh = fD.mesh();

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

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

            vectorField Md = mesh.C().boundaryField()[patchi] - CofR_;

            scalarField sA = mag(Sfb[patchi]);

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

            fm.first().first() += sum(fN);
            fm.second().first() += sum(Md ^ fN);

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

            fm.first().second() += sum(fT);
            fm.second().second() += sum(Md ^ fT);
        }
    }
    else
    {
        const volVectorField& U = obr_.lookupObject<volVectorField>(UName_);
        const volScalarField& p = obr_.lookupObject<volScalarField>(pName_);

        const fvMesh& mesh = U.mesh();

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

        tmp<volSymmTensorField> tdevRhoReff = devRhoReff();
        const volSymmTensorField::GeometricBoundaryField& devRhoReffb
            = tdevRhoReff().boundaryField();

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

            vectorField Md = mesh.C().boundaryField()[patchi] - CofR_;

            vectorField pf = Sfb[patchi]*p.boundaryField()[patchi];

            fm.first().first() += rho(p)*sum(pf);
            fm.second().first() += rho(p)*sum(Md ^ pf);

            vectorField vf = Sfb[patchi] & devRhoReffb[patchi];

            fm.first().second() += sum(vf);
            fm.second().second() += sum(Md ^ vf);
        }
    }

    reduce(fm, sumOp());

    return fm;
}

Can anyone by this understand, how OpenFOAM calculates these forces?

Thanks a lot for all your help!

Hi, for calculating total Forces in OpenFOAM, integrated the pressure Field over the surfaces and projection of this forces in drag or lift direction Division Pdyn*Aref gives the cd and cl Coeff

_____
Rasoul

[malhar]

hello every1,
i am also new to openfoam. i have done similar sim ulation for flow acrodss a cylinder. i get drag n lift forces correctly, i.e accor ding to vortex shedding i am getting variation in lift forces. but the coeff of drag n lift Cd n Cl, i am getting them constant throughout. i am unable ti digest this contrasting behavour.
please guide me thr this.

thanks n regards malhar.

[malaboss]

Originally Posted by desert_1250

Hi, for calculating total Forces in OpenFOAM, integrated the pressure Field over the surfaces and projection of this forces in drag or lift direction Division Pdyn*Aref gives the cd and cl Coeff

Hi,
by saying this, do you mean that the pressure value used is the one on the boundary of the cells which touch the surface on which drag and lift force are applied ?

I am a new Foamer but i think that we forget the viscous force if OpenFOAM only does this. In the code, there is a tangential force calculated, but i cannot say if it concerns the viscous force or not. Does anyone have an idea ?

[Tushar@cfd]

Understanding of the terms...

I am actually looking out for the term mentioned below in the previous posts on this thread:

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

Where is the term located in the case directory?

Correct me if I am wrong

Thanks in Advance!

[Tushar@cfd]

Oh, it's clear the term is for compressible flows.

Please, excuse for previous post. I did mistake it for incompressible flows.

[ck7]

hi,

https://www.youtube.com/watch?time_c...ture=emb_title

The above link provides a way to calculate drag and lift using ParaView.