How can I calculate the

- mass flux averaged pressure for a boundary region
- mass flow rate for a region (sum, inflow, outflow)

during the calculation?

Have a look at this code example:

/********************************************/
scalar volumeFluxWeightedPressureInlet = 0.0;
scalar volumeFluxInlet = 0.0;
word inletPatchName = "inlet";

//access boundary elements
const fvPatchList& patches = mesh.boundary();

//loop over boundaries
forAll(patches, patchI)
{
<BLOCKQUOTE> const fvPatch cPatch = patches[patchI];

//check boundary name
if(cPatch.name() == inletPatchName)
{
<BLOCKQUOTE> //Access boundary face area vectors
vectorField inletAreaVectors = cPatch.faceAreas();

//loop over all faces of selected
//boundary
forAll(cPatch, faceI)
{
<BLOCKQUOTE> //calculate boundary face flux
scalar cFaceFlux
= U.boundaryField()[patchI][faceI] & inletAreaVectors[faceI];

//sum face fluxes
volumeFluxInlet += cFaceFlux;

//sum flux weighted face pressures
volumeFluxWeightedPressureInlet += p.boundaryField()[patchI][faceI] * cFaceFlux;</BLOCKQUOTE></BLOCKQUOTE><BLOCKQUOTE> }</BLOCKQUOTE><BLOCKQUOTE> volumeFluxWeightedPressureInlet \= volumeFluxInlet;</BLOCKQUOTE></BLOCKQUOTE><BLOCKQUOTE> }</BLOCKQUOTE>}
/**********************************************/


Here "mesh" is the fvMesh grid, U is the velocity field (vector), p is the pressure field and "&" is the dot product operator.

The same general procedure can be used to access, modify and otherwise manipulate all types of boundary field values. For example, it can also be used to specify inlet velocities through face centre locations and or time dependent functions.

See "src/OpenFOAM/lnInclude/fvPatch.H" and its base classes in the OpenFOAM directory or the Doxygen documentation for more relevant access functions.

Eugene

Hi Eugene,

I tried to get it running but I still get one error message:

error: 'const class Foam::fvPatch' has no member named 'faceAreas'

I think that I got an idea what you are doing in the example, but I'm a total c++ beginner ;-) I will also send you the example via email.

Thanks very much for your help

Jörn

Try the function Sf() instead (faceAreas is a polyPatch member function, which I mistakenly thought would be inherited by fvPatch).

For a C++ beginner, the Doxygen documentation can be extremely valuable; providing an easy to understand summary of class member functions and inheritance. It also has a built-in search engine, which saved me loads of time when I first started to use Foam and was still struggling with the basics of C++.

In response to your e-mail, yes "rho.boundaryField()[patchI][faceI]" will give you the density at a specified face on the specified boundary patch.

Eugene

I don't agree with Eugene's example as a method of calculating the patch mass flow rate. It is the flux field phi which contains the mass flux for compressible codes or the volumetic flux for incompressible codes. Also it is not usually neccessary to loop over faces; there is a large set of field functions included in FOAM which sould be used in preference to low-level looping where ever possible. Thus the inlet mass flux can simply be obtained from:

scalar massFlux = sum(phi.boundaryField()[inletPatchi]);

where the inletPatchi is the index of the inlet patch which can be obtined by searching as in Eugene's example or from

label inletPatchi = mesh.boundaryMesh().findPatchID(nameOfInlet);

where nameOfInlet is the name of the inlet patch which could be provided from a dictionary etc.

Thanks Eugene and Henry,

after removing some typos and changing to the Sf() function the compilation and calculation work perfectly well.

I already had a look into the Doxygen documentation but for a c++ beginner it looks a bit like "chinese backwards" ;-)

For those who want to use it I will include the running version here:

/********************************************/

scalar massFluxWeightedPressureInlet = 0.0 ;
scalar massFluxInlet = 0.0 ;
word inletPatchName = "INLET";

//access boundary elements
const fvPatchList& patches = mesh.boundary();

//loop over boundaries
forAll(patches, patchI)
{

const fvPatch& cPatch = patches[patchI];

//check boundary name
if(cPatch.name() == inletPatchName)
{

//Access boundary face area vectors
const vectorField& inletAreaVectors = cPatch.Sf();

//loop over all faces of selected
//boundary
forAll(cPatch, faceI)
{

//calculate boundary face mass flux
scalar cFaceFlux
= U.boundaryField()[patchI][faceI] & inletAreaVectors[faceI] *
rho.boundaryField()[patchI][faceI];

//sum face mass fluxes
massFluxInlet += cFaceFlux;

//sum mass flux weighted face pressures
massFluxWeightedPressureInlet += p.boundaryField()[patchI][faceI] * cFaceFlux;

}

massFluxWeightedPressureInlet /= massFluxInlet;

}

}

Info<< "MassFlowRate: " << massFluxInlet << " kg/s" << endl;
Info<< "Pressure : " << massFluxWeightedPressureInlet << " Pa " << endl;

/**********************************************/

I still think this is an absolutely terrible way of doing things in FOAM. I posted a much more elegant and appropriate way of calculating the patch-flux and I must warn you again that for outlets in particular the flux obtained from the velocity in this way does not obey continuity.

Ok Henry,

I replaced the

//calculate boundary face mass flux
scalar cFaceFlux
= U.boundaryField()[patchI][faceI] & inletAreaVectors[faceI] *
rho.boundaryField()[patchI][faceI];

by

scalar cFaceFlux
= phi.boundaryField()[patchI][faceI];

Your suggestion works for calculating the flux without looping the faces. But is there a way to calculate the mass flow averaged pressure without this loop?

PS: the pressure difference in now below 2% with ke and ud.

You can do nearly anything in FOAM without looping. Imagine you have a field psi and want a flux-weighted average of it over patch patchi:

sum(phi.boundaryField()[patchi]*psi.boundaryField()[patchi])/sum(phi.boundaryField()[patchi])

I guess that "...nearly anything in FOAM without looping" is supposed to be read as "...nearly anything in FOAM without explicit looping"

I am using the following code to summarise the flux at each boundary patch, along with the net, after each timestep.

{
// loop through all boundary patches
int patchItr;
label patchIndex;
scalar flux, netFlux = 0.0;
forAll (mesh.boundaryMesh(), patchItr)
{

// get the patch index
patchIndex = mesh.boundaryMesh()[patchItr].index();

// compute the mass flux and net flux
flux = sum(phi.boundaryField()[patchIndex]);
netFlux += flux;

// compute and display the mass flux with patch name
Info<< "Mass flux at "
<< mesh.boundaryMesh()[patchIndex].name()
<< " = " << flux
<< endl;
}
Info<< "Net mass flux = " << netFlux << endl << endl;
}

This works fine, but in parallel it can't resolve some of the boundary names that are not available on the master processor (giving a name of "procBoundary0to1" for example). How can the master processor get access to these names via the mesh.boundaryMesh()[patchIndex].name() statement? (I couldn't find a solution in the Doxygen notes). However, the flux appears to be calculated fine.

"procBoundary0to1" is the processor-processor "boundary" patch which isn't a true boundary, it's the internal patch created by the decomposition and probably not important to your flux calculation. In order to sum over the "real" boundary patches in parallel you will need to test that the patches you are interested are on the processor and after the loop do a reduce operation to sum up the contributions from the processors to the boundary regions. Have a look through src/cfdTools/adjustPhi/adjustPhi.C which uses the logic you need to sum inlet and outlet fluxes to correct mass-balance for steady-state cases.

Thanks for that Henry, I have made those modifications and am using the reduce command to sum the contributions from each processor. Is there a way to get a list of all boundaries on the original mesh while running in parallel? Or is this not possible due to running on decomposed grids in seperate processes? At the moment I am defining a patch list explicity, but I could potentially parse the original constant/polyMesh/boundary file to automate this (assuming the files are available during runtime).

You can send the information about the boundary patches to the master, sort the lists and send them back. However, the patch numbering on the processors will be different because each only hold the list of it's own patches so I am not sure how this complete list will help you.

[sinclair]

I think I've got a solution that works in both serial and parallel, and uses the communcation and sorting you suggested.

Firstly, you need to #include "SortableList.H" as one of the top level includes in a solver application.

Then #include "buildGlobalBoundaryList.H" once "createMesh.H" and other initialisations have been run in the main function. (This only needs to be run once.) The file is attached below.

buildGlobalBoundaryList.H

The result of this is a wordList variable called "globalPatchNames" that is available on each processor, and contains a sorted list of all real patches in the model.

To compute mass fluxes after each timestep, #include "computeMassFlux.H" (see attachment).

computeMassFlux.H

I've tested this in both serial and parallel, so hopefully it works in other people's cases.

[henry]

Thanks for these pieces of sample code, it all seems fine except for one point, did you need to include Pstream.H? If so it should be included outside the "main" function at the top of the code although it is usually already included via some other includes.

[sinclair]

You are right about Pstream.H Henry, I've made some changes to the above code. buildGlobalBoundaryList.H is now a function that returns a wordList of patch names rather than inline code, so to avoid any variable redefinition problems. This also makes the headers a bit cleaner. However, when I perform a sort on a SortableList with OpenFOAM 1.1, I get a segmentation fault (this worked in 1.0.2). Do you know why this is so?

[henry]

The difference between 1.0.2 and 1.1 with respect to sortable list is that I have changed all calls to the C-library qsort and replaced them the equivalent STL sort calls which is certainly a more elegant method but we have had problems with it in the past.

Could you please send us the code that produces the failure and we will look into it.

[sinclair]

Sure, the new version of "buildGlobalBoundaryList.H" is attached below:

buildGlobalBoundaryList.H

This should be included outside the main function after top level includes like "fvCFD.H". It can then be called inside the main function like so,

wordList globalBoundaryList;
globalBoundaryList = buildGlobalBoundaryList(mesh);

Then, mass fluxes can be calculated via the following include (similar to the previous example).

computeMassFlux.H

I've commented out the sort() operation in "buildGlobalBoundaryList.H" for now.

[mattijs]

Hello Jarrod,

we located found the problem in SortableList. It is because setSize when it grows the list does not adapt the list of indices. Until we fix it you can work around it by creating your list big enough and shrinking it (i.e. setSize with something smaller than the current size)

Mattijs