CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Post-Processing (http://www.cfd-online.com/Forums/openfoam-post-processing/)
-   -   Power draw of a mixer (http://www.cfd-online.com/Forums/openfoam-post-processing/74659-power-draw-mixer.html)

erklaerbaer April 6, 2010 04:33

Power draw of a mixer
 
Hi all,

my goal is to evaluate the power draw of a mixer (mixer3D) via the integration of the forces acting on the moving or non-moving parts over the relevant patch.

I am an absolute Newbie in OpenFOAM and threw together some Code from the wallShearStree, wallGradU, patchIntegrate utilities and some things i read in the forum.

I try to get the velocity gradient with

Code:

IOobject Uheader
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ
        );
 
Info<< "    Reading U" << endl;
volVectorField U(Uheader, mesh);
Info<< "    Calculating wallGradU" << endl;
volVectorField wallGradU
            (   
        IOobject
                (
                    "wallGradU",
                    runTime.timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::AUTO_WRITE
                ),
        mesh,
              dimensionedVector
                (
                    "wallGradU",
                    U.dimensions()/dimLength,
                    vector::zero
                )
            );
       
    forAll(wallGradU.boundaryField(), patchi)
            {
                wallGradU.boundaryField()[patchi] =
                    U.boundaryField()[patchi].snGrad();
   
            }

and to calculate the sheer rate with
Code:

Info<< "    Calculating Stress" << endl;
        volScalarField wallForces
        (
            IOobject
            (
                "wallForces",
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            //dynVis*sqrt(sqr(wallGradU.component(vector::X)) + sqr(wallGradU.component(vector::Y))+sqr(wallGradU.component(vector::Z)))
        dynVis*mag(wallGradU)
        );

      wallForces.write();

If i am correct the value for tau is at this point calculated and stored in wallForces. Am i right?

Now i want to integrate tau over the faces of the patch, or not?

Code:

IOobject fieldHeader
        (
            "wallForces",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ
        );

        // Check field exists
        if (fieldHeader.headerOk())
        {
            mesh.readUpdate();

            label patchi = mesh.boundaryMesh().findPatchID(patchName);
            if (patchi < 0)
            {
                FatalError
                    << "Unable to find patch " << patchName << nl
                    << exit(FatalError);
            }

const polyPatch& cPatch = mesh.boundaryMesh()[patchi];
const surfaceScalarField& magSf = mesh.magSf();
   
    scalar Power = 0;
    volScalarField field(fieldHeader, mesh);

    forAll(cPatch, facei)
    {
        Power += magSf.boundaryField()[patchi][facei] * field.boundaryField()[patchi][facei];
    }
    Info<< "    Power = " << Power << endl;

Now the biggest Problem is that the calculation is not finished because the Integration of tau over the patch area gives me a force but not Energy or even Power.
Can someone give me a hint how i bring this "element" into my code? Or point me to a source where i can find the formula to calcualte the power draw via this method.

Thanks in advance.
Andi

holger_marschall April 6, 2010 14:51

Hi Andi,

welcome to this board...

You might want to have a look at
Handbook of Industrial Mixing (page 314)
Published Online: 30 Jan 2004
http://www3.interscience.wiley.com/c...home/107583875
Editor(s): Edward L. Paul, Victor A. Atiemo-Obeng, Suzanne M. Kresta
Print ISBN: 9780471269199 Online ISBN: 9780471451457
DOI: 10.1002/0471451452

"The power delivered to the fluid is the product of the impeller speed, 2πN, in rad/s, and torque, τ, which is obtained by integration of the pressure on the impeller blade: P = 2πNτ"

and/or
CFD Simulation of Flow in a Baffled Stirred Tank
GL Lane, PTL Koh
Chemeca 97, 25th Australian and New Zealand Chemical Engineering Conference, Rotorua, NZ, September/October 1997, Paper FL2B, 9 pp
http://www.cfd.com.au/cfd_conf97/papers/lan035.pdf
This should give a first approximate number (if one neglects the shear). Hope this helps.

best,

erklaerbaer April 7, 2010 15:56

Thank you .

Unfortunately I am stuck at the next phase. But now it is more a problem of the OpenFOAM Programming syntax i dont really geht.

I try to get torque by summing up the
Area of every Face * Pressure on that face * Distance to the Center (= lever arm)

My Syntax for that is

Code:

##Read the calculated pressure from the corresponding file
IOobject pressureHeader
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ
        );
volScalarField pressure(pressureHeader, mesh);

#Find IF for the Patch, in that case the Impeller
label patchi = mesh.boundaryMesh().findPatchID(patchName);
const polyPatch& cPatch = mesh.boundaryMesh()[patchi];

#I guess this reads the areasizes
const surfaceScalarField& magSf = mesh.magSf();


forAll(cPatch, facei)
    {
        Torque +=
magSf.boundaryField()[patchi][facei] *
mag(pressure.boundaryField()[patchi][facei])*
Foam::sqrt(
    sqr(pressure.mesh().boundary()[patchi].Cf()[facei][0])+
    sqr(pressure.mesh().boundary()[patchi].Cf()[facei][1])
      );
}

magSf.boundaryField()[patchi][facei] = area of the face?
mag(pressure.boundaryField()[patchi][facei]) = pressure on that face?
pressure.mesh().boundary()[patchi].Cf()[facei][0] = The x-Coordinate of the Center of the face?

The results i get are to low by arround 1E-2 so i guess i have a missunderstanding in the Syntax of OpenFOAM, if someone could lighten me up i would be very grateful.

Andi

erklaerbaer April 8, 2010 05:38

1 Attachment(s)
If somebody is interested:

I wrote a little Tool for turbulent flows which calculates the powerdraw from the Epsilon-Term,
P = rho * Integral(Epsilon dV)

andi

maximeg June 14, 2010 15:43

Hi!

I'm a beginner at OpenFOAM and I'd like some help if possible.

I have to calculate the circulation in a vortex. My code is as follows:

Code:

#include "calc.H"
#include "fvc.H"

void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
{

    volScalarField magVorticity        // Imports the results of the vorticity of the vortex
    (                   
        IOobject           
        (               
            "magVorticity",       
            runTime.timeName(),       
            mesh,           
            IOobject::MUST_READ,   
            IOobject::AUTO_WRITE   
        ),               
        mesh               
    );     

    scalar circulation = 0.0;
    scalar r = 0.0;
    const volScalarField& magSf = mesh.magSf();

    forAll(magVorticity, cellI)
    {

        if (distance < r[cellI])        // If the radius (distance) stays lower than the critical radius (r)...
        {
            circulation += magSf.magVorticity[cellI];    // ...calculation of the circulation.
        }

    }

}

The error message I get in the terminal is as follows:

Code:

SOURCE=MaxVorCirc.C ;  g++ -m64 -Dlinux64 -DWM_DP -Wall -Wno-strict-aliasing -Wextra -Wno-unused-parameter -Wold-style-cast -Wnon-virtual-dtor -O3  -DNoRepository -ftemplate-depth-40 -I/usr/local/OpenFOAM/OpenFOAM-1.6/src/postProcessing/postCalc -I/usr/local/OpenFOAM/OpenFOAM-1.6/src/finiteVolume/lnInclude -IlnInclude -I. -I/usr/local/OpenFOAM/OpenFOAM-1.6/src/OpenFOAM/lnInclude -I/usr/local/OpenFOAM/OpenFOAM-1.6/src/OSspecific/POSIX/lnInclude  -fPIC -c $SOURCE -o Make/linux64GccDPOpt/MaxVorCirc.o
MaxVorCirc.C: In function ‘void Foam::calc(const Foam::argList&, const Foam::Time&, const Foam::fvMesh&)’:
MaxVorCirc.C:89: error: invalid types ‘Foam::scalar[Foam::label]’ for array subscript
MaxVorCirc.C:94: error: ‘magSf’ was not declared in this scope
make: *** [Make/linux64GccDPOpt/MaxVorCirc.o] Error 1

Any help would be very appreciated.

Thank you in advance!

Maxime


All times are GMT -4. The time now is 12:27.