CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Running, Solving & CFD

Calculation of integral boundary values

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

Like Tree2Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   January 12, 2005, 04:01
Default How can I calculate the -
  #1
Joern Beilke (Beilke)
Guest
 
Posts: n/a
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?
  Reply With Quote

Old   January 12, 2005, 06:54
Default Have a look at this code examp
  #2
Eugene de Villiers (Eugene)
Guest
 
Posts: n/a
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
  Reply With Quote

Old   January 12, 2005, 15:32
Default Hi Eugene, I tried to get
  #3
Joern Beilke (Beilke)
Guest
 
Posts: n/a
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
  Reply With Quote

Old   January 12, 2005, 19:22
Default Try the function Sf() instead
  #4
Eugene de Villiers (Eugene)
Guest
 
Posts: n/a
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
  Reply With Quote

Old   January 13, 2005, 04:50
Default I don't agree with Eugene's e
  #5
Henry Weller (Henry)
Guest
 
Posts: n/a
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.
sharonyue likes this.
  Reply With Quote

Old   January 13, 2005, 08:33
Default Thanks Eugene and Henry, a
  #6
Joern Beilke (Beilke)
Guest
 
Posts: n/a
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;

/**********************************************/
  Reply With Quote

Old   January 13, 2005, 08:46
Default I still think this is an abso
  #7
Henry Weller (Henry)
Guest
 
Posts: n/a
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.
  Reply With Quote

Old   January 13, 2005, 09:28
Default Ok Henry, I replaced the
  #8
Joern Beilke (Beilke)
Guest
 
Posts: n/a
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.
  Reply With Quote

Old   January 13, 2005, 10:15
Default You can do nearly anything in
  #9
Henry Weller (Henry)
Guest
 
Posts: n/a
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])
  Reply With Quote

Old   January 13, 2005, 11:08
Default I guess that "...nearly anyth
  #10
Nils Smeds (Smeds)
Guest
 
Posts: n/a
I guess that "...nearly anything in FOAM without looping" is supposed to be read as "...nearly anything in FOAM without explicit looping"
  Reply With Quote

Old   March 3, 2005, 23:46
Default I am using the following code
  #11
Jarrod Sinclair (Sinclair)
Guest
 
Posts: n/a
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.
  Reply With Quote

Old   March 4, 2005, 03:39
Default "procBoundary0to1" is the pro
  #12
Henry Weller (Henry)
Guest
 
Posts: n/a
"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.
  Reply With Quote

Old   March 7, 2005, 03:11
Default Thanks for that Henry, I have
  #13
Jarrod Sinclair (Sinclair)
Guest
 
Posts: n/a
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).
  Reply With Quote

Old   March 7, 2005, 04:48
Default You can send the information
  #14
Henry Weller (Henry)
Guest
 
Posts: n/a
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.
  Reply With Quote

Old   March 8, 2005, 22:02
Default I think I've got a solution th
  #15
New Member
 
Jarrod Sinclair
Join Date: Mar 2009
Posts: 6
Rep Power: 8
sinclair is on a distinguished road
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.
mm.abdollahzadeh likes this.
sinclair is offline   Reply With Quote

Old   March 9, 2005, 04:51
Default Thanks for these pieces of sam
  #16
Senior Member
 
Join Date: Mar 2009
Posts: 854
Rep Power: 13
henry is on a distinguished road
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.
henry is offline   Reply With Quote

Old   March 16, 2005, 20:07
Default You are right about Pstream.H
  #17
New Member
 
Jarrod Sinclair
Join Date: Mar 2009
Posts: 6
Rep Power: 8
sinclair is on a distinguished road
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?
sinclair is offline   Reply With Quote

Old   March 17, 2005, 04:31
Default The difference between 1.0.2 a
  #18
Senior Member
 
Join Date: Mar 2009
Posts: 854
Rep Power: 13
henry is on a distinguished road
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.
henry is offline   Reply With Quote

Old   March 17, 2005, 09:08
Default Sure, the new version of "buil
  #19
New Member
 
Jarrod Sinclair
Join Date: Mar 2009
Posts: 6
Rep Power: 8
sinclair is on a distinguished road
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.
sinclair is offline   Reply With Quote

Old   March 17, 2005, 09:58
Default Hello Jarrod, we located fo
  #20
Super Moderator
 
Mattijs Janssens
Join Date: Mar 2009
Posts: 1,416
Rep Power: 16
mattijs is on a distinguished road
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
mattijs 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
Surface integral calculation by undefined area Bloshchitsyn Vladimir CFX 0 November 18, 2007 23:19
Evaluation of boundary integral aunola Main CFD Forum 0 May 8, 2006 07:54
Boundary Integral Method CFDtoy Main CFD Forum 0 October 5, 2004 13:08
Boundary integral parameters StarCD CD-adapco 2 May 24, 2004 14:33
Calculation of lagrangian integral time scale R.Sripriya FLUENT 4 June 12, 2003 12:49


All times are GMT -4. The time now is 05:41.