CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Mass weighted average

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   July 31, 2013, 12:46
Default Mass weighted average
  #1
Member
 
Patrick Wollny
Join Date: Apr 2010
Posts: 58
Rep Power: 7
Pat84 is on a distinguished road
Hiho,

Im modifying the tool patchAverage with the intention to calculate the mass weighted average, but I get an error while compiling the code. The reason is the "&" operator, which should make a inner product of the patch normal vector and the velocity vector at the boundary.

Here is the code (the line with the vector operation is marked in red):



Code:
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright held by original author
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software; you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the
    Free Software Foundation; either version 2 of the License, or (at your
    option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM; if not, write to the Free Software Foundation,
    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

Application
    massAverage

Description
    Calculates the average of the specified field over the specified patch.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:

int main(int argc, char *argv[])
{
    timeSelector::addOptions();
    argList::validArgs.append("fieldName");
    argList::validArgs.append("patchName");
#   include "setRootCase.H"
#   include "createTime.H"
    instantList timeDirs = timeSelector::select0(runTime, args);
#   include "createMesh.H"

    word fieldName(args.additionalArgs()[0]);
    word patchName(args.additionalArgs()[1]);
    //dimensionedScalar def_rho("def_rho",dimensionSet(1,-3,0,0,0,0,0),args.additionalArgs()[2]);
    //dimensionedScalar def_rho("def_rho",dimensionSet(1,-3,0,0,0,0,0),1.226);

    forAll(timeDirs, timeI)
    {
        runTime.setTime(timeDirs[timeI], timeI);
        Info<< "Time = " << runTime.timeName() << endl;

        IOobject fieldHeader
        (
            fieldName,
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ
        );
    
    IOobject U_header
        (
            "U",
            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);
            }

            if (fieldHeader.headerClassName() == "volScalarField")
            {
                Info<< "    Reading volScalarField " << fieldName << endl;
                volScalarField field(fieldHeader, mesh);
        volVectorField U(U_header, mesh);
        volScalarField UA = U.boundaryField()[patchi].patchInternalField() & mesh.Sf().boundaryField()[patchi];
                //volScalarField UA = U.boundaryField()[patchi] & mesh.Sf().boundaryField()[patchi];
        
                scalar total_mass = gSum(UA.boundaryField()[patchi] * def_rho);
                scalar sumField = 0;

                if (total_mass > 0)
                {
                    sumField = gSum
                    (
                        UA.boundaryField()[patchi] * def_rho
                      * field.boundaryField()[patchi]
                    ) / total_mass;
                }

                Info<< "    Average of " << fieldName << " over patch "
                    << patchName << '[' << patchi << ']' << " = "
                    << sumField << endl;

                Info<< "    mesh.Sf().boundaryField()[patchi]   " << mesh.Sf().boundaryField()[patchi] << "\n"
                    << "    U.boundaryField()[patchi]           " << U.boundaryField()[patchi] << endl; 


            }
            else
            {
                FatalError
                    << "Only possible to average volScalarFields "
                    << nl << exit(FatalError);
            }
        }
        else
        {
            Info<< "    No field " << fieldName << endl;
        }

        Info<< endl;
    }

    Info<< "End\n" << endl;

    return 0;
}

// ************************************************************************* //
volVecorFields are more lists than vectors. But why does the "*" operator work and how can I perform vector and tensor operations?

Any ideas?
Pat84 is offline   Reply With Quote

Old   July 31, 2013, 17:51
Default
  #2
Member
 
Patrick Wollny
Join Date: Apr 2010
Posts: 58
Rep Power: 7
Pat84 is on a distinguished road
The code is now running. It was a silly mistake.

This is the code:

Code:
Application
    massAverage

Description
    Calculates the average of the specified field over the specified density.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:

int main(int argc, char *argv[])
{
    timeSelector::addOptions();
    argList::validArgs.append("fieldName");
    argList::validArgs.append("patchName");
#   include "setRootCase.H"
#   include "createTime.H"
    instantList timeDirs = timeSelector::select0(runTime, args);
#   include "createMesh.H"

    word fieldName(args.additionalArgs()[0]);
    word patchName(args.additionalArgs()[1]);
    //dimensionedScalar def_rho("def_rho",dimensionSet(1,-3,0,0,0,0,0),args.additionalArgs()[2]);
    //dimensionedScalar def_rho("def_rho",dimensionSet(1,-3,0,0,0,0,0),scalar(1.226));
    
    volScalarField def_rho
        (
          IOobject
            (
              "def_rho",
              runTime.timeName(),
              mesh,
              IOobject::NO_READ,
              IOobject::NO_WRITE
            ),
            mesh,
            dimensionedScalar
            (
              "def_rho",
               dimensionSet(1,-3,0,0,0,0,0),
                       scalar(1.226)
                    )
                 );

    forAll(timeDirs, timeI)
    {
        runTime.setTime(timeDirs[timeI], timeI);
        Info<< "Time = " << runTime.timeName() << endl;

        IOobject fieldHeader
        (
            fieldName,
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ
        );
    
    IOobject U_header
        (
            "U",
            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);
            }

            if (fieldHeader.headerClassName() == "volScalarField")
            {
                Info<< "    Reading volScalarField " << fieldName << endl;
                volScalarField field(fieldHeader, mesh);
        volVectorField U(U_header, mesh);
        
        volScalarField UA
        (
          IOobject
            (
              "UA",
              runTime.timeName(),
              mesh,
              IOobject::NO_READ,
              IOobject::NO_WRITE
            ),
            mesh,
            dimensionedScalar
            (
              "UA",
               dimensionSet(0,2,-1,0,0,0,0),
                       0
                    )
                 );
        
        
        
        UA.boundaryField()[patchi] = (U.boundaryField()[patchi] & mesh.Sf().boundaryField()[patchi]);
        
                scalar total_mass = gSum(UA.boundaryField()[patchi] * def_rho.boundaryField()[patchi]);
                scalar sumField = 0;

                if (total_mass > 0)
                {
                    sumField = gSum
                    (
                        UA.boundaryField()[patchi] * def_rho
                      * field.boundaryField()[patchi]
                    ) / total_mass;
                }

                Info<< "    Average of " << fieldName << " over patch "
                    << patchName << '[' << patchi << ']' << " = "
                    << sumField << endl;

            }
            else
            {
                FatalError
                    << "Only possible to average volScalarFields "
                    << nl << exit(FatalError);
            }
        }
        else
        {
            Info<< "    No field " << fieldName << endl;
        }

        Info<< endl;
    }

    Info<< "End\n" << endl;

    return 0;
}

// ************************************************************************* //
The red lines do not satisfying me, because I need to allocate a whole list with the same scalar for multiplication with the scalar UA. Is it possible to allocate only one variable def_rho ( like in the out-comment line with dimensonedScalar) for this multiplication?

Best regards,
Pat
Pat84 is offline   Reply With Quote

Old   August 1, 2013, 06:46
Default
  #3
Member
 
Patrick Wollny
Join Date: Apr 2010
Posts: 58
Rep Power: 7
Pat84 is on a distinguished road
Since sum(scalar*rho*U*A)/sum(rho*U*A) = sum(scalar*U*A)/sum(U*A) for rho = const
and the inner product of U and A is saved in phi, I have changed my code:

Code:
int main(int argc, char *argv[])
{
    timeSelector::addOptions();
    argList::validArgs.append("fieldName");
    argList::validArgs.append("patchName");
#   include "setRootCase.H"
#   include "createTime.H"
    instantList timeDirs = timeSelector::select0(runTime, args);
#   include "createMesh.H"

    word fieldName(args.additionalArgs()[0]);
    word patchName(args.additionalArgs()[1]);

    forAll(timeDirs, timeI)
    {
        if(timeI == 0) continue;
      
        runTime.setTime(timeDirs[timeI], timeI);
        Info<< "Time = " << runTime.timeName() << endl;

        IOobject fieldHeader
        (
            fieldName,
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ
        );
    
    IOobject phi_header
        (
            "phi",
            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);
            }

            if (fieldHeader.headerClassName() == "volScalarField")
            {
                Info<< "    Reading volScalarField " << fieldName << endl;
                volScalarField field(fieldHeader, mesh);
        surfaceScalarField phi(phi_header, mesh);

        
                scalar total_mass = gSum(phi.boundaryField()[patchi]);
                scalar sumField = 0;

                if (total_mass > 0)
                {
                    sumField = gSum
                    (
                        phi.boundaryField()[patchi]
                      * field.boundaryField()[patchi]
                    ) / total_mass;
                }

                Info<< "    Average of " << fieldName << " over patch "
                    << patchName << '[' << patchi << ']' << " = "
                    << sumField << endl;

            }
            else
            {
                FatalError
                    << "Only possible to average volScalarFields "
                    << nl << exit(FatalError);
            }
        }
        else
        {
            Info<< "    No field " << fieldName << endl;
        }

        Info<< endl;
    }

    Info<< "End\n" << endl;

    return 0;
}
It also is now more accurate since it uses the fluxes stored in the face centers.

You can now use it for openFoam1.6 and openFoam1.6-ext.
Attached Files
File Type: gz massAverage.tar.gz (1.7 KB, 17 views)

Last edited by Pat84; August 5, 2013 at 10:36.
Pat84 is offline   Reply With Quote

Old   April 16, 2014, 04:20
Default
  #4
Senior Member
 
romant's Avatar
 
Roman Thiele
Join Date: Aug 2009
Location: Stockholm, Sweden
Posts: 359
Rep Power: 11
romant is on a distinguished road
Hej,

I think the tool is ready needed. However, I was wondering if you can update it to be compiled in OF 2.3? I just checked to compile it. I get only one error message when compiling and I commented out one line in order to get it to compile. Line 53, what is it for?
Code:
        // if(timeDirs[timeI] == 0) continue;
I then wanted to apply it to one of my cases. It seems that your tool runs through all patches and not just the one specified, because it checks the wall patches that I didn't specify in the call for massAverage T outlet as well.

It only "knows" the standard patches, and does not know about the turbulence patch types. I think you can change that by actually including them in the compiling process. For example in my case it does not recognize the patch type compressible::turbulentHeatFluxTemperature.
__________________
~roman
romant is offline   Reply With Quote

Old   July 22, 2014, 15:13
Default
  #5
New Member
 
Jens
Join Date: Oct 2013
Posts: 2
Rep Power: 0
Jens123 is on a distinguished road
Hi,

good Tool. Ive searched exactly for this. I am very lucky to find it

In my opinion Line 53 should only skip the Time "Zero" and jump directly to the next time step. Because in Time "Zero" there is no "phi" file and nobody need the information about this time step.
Jens123 is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Area weighted or mass weighted average SAM Main CFD Forum 24 July 8, 2015 06:15
y+ and u+ values with low-Re RANS turbulence models: utility + testcase florian_krause OpenFOAM 108 June 9, 2015 08:13
mass weighted average temperature along a curve preetam69 Visualization & Post-Processing 1 January 14, 2014 12:33
mass weighted average temperature along a curve preetam69 FLUENT 1 July 5, 2013 02:48
Diffference between mass weighted average and Area giogio FLUENT 4 September 25, 2012 01:36


All times are GMT -4. The time now is 00:10.