CFD Online Logo CFD Online URL
Home > Forums > OpenFOAM Post-Processing

Volume averaging in cell

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

LinkBack Thread Tools Display Modes
Old   February 17, 2014, 14:44
Post Volume averaging in cell
Vincent Leroy
Join Date: Jul 2012
Location: Bordeaux, France
Posts: 38
Rep Power: 5
leroyv is on a distinguished road
Dear OpenFOAM users,

I solved a flow problem over and array of cylinders:

I would like to get the volume average of the velocity and pressure fields over a limited, moving volume (typically the size of a single cell, see the mask on the screenshot). The resulting field would be continuously defined over the entire domain, including in the cylinders. This is field volume averaging / spatial frequency filtering, in the sense of Whitaker (1999) and is somehow like the spatial averaging operation used for LES theory.

Do you think it is possible to do this while solving the problem or as a post-processing operation? I didn't find a way to apply such a moving filter in Paraview or in the sample utility, and I don't know how to program this in a solver.

Last edited by leroyv; February 18, 2014 at 13:25. Reason: added more precision
leroyv is offline   Reply With Quote

Old   February 18, 2014, 12:57
Post Update
Vincent Leroy
Join Date: Jul 2012
Location: Bordeaux, France
Posts: 38
Rep Power: 5
leroyv is on a distinguished road
Additional (maybe simpler) question: How to export the fields to that kind format:
x y z Ux Uy Uz p
with (x,y,z) coordinates on a cartesian grid. Just doing that would be sufficient, the computation of the average being possible as a post-processing operation.

Last edited by leroyv; February 21, 2014 at 10:04. Reason: Removed quote
leroyv is offline   Reply With Quote

Old   February 27, 2014, 11:09
Smile Solution (non optimal)
Vincent Leroy
Join Date: Jul 2012
Location: Bordeaux, France
Posts: 38
Rep Power: 5
leroyv is on a distinguished road
After a week of hacking (or so), I put together a solution. I forked simpleFoam and added this piece of code in it:

    Info<<"Computing average velocity field\n" << endl;

    // Constants
    IOdictionary filterProperties(

    // This is AWFUL hacking
    scalar  filterRadius = dimensionedScalar(filterProperties.lookup("filterRadius")).value(),
            filterXMin   = dimensionedScalar(filterProperties.lookup("filterXMin")).value(),
            filterXMax   = dimensionedScalar(filterProperties.lookup("filterXMax")).value(),
            filterYMin   = dimensionedScalar(filterProperties.lookup("filterYMin")).value(),
            filterYMax   = dimensionedScalar(filterProperties.lookup("filterYMax")).value();
    int     nX           = dimensionedScalar(filterProperties.lookup("nX")).value(),
            nY           = dimensionedScalar(filterProperties.lookup("nY")).value();

    scalar  hX = (filterXMax - filterXMin) / int(nX-1),
            hY = (filterYMax - filterYMin) / int(nY-1);
    const volVectorField & cellCenterField = mesh.C();
    const DimensionedField<scalar, volMesh> & cellVolumeField = mesh.V();

    // Variable initialization
    scalar normInfty = 0.;
    scalar filterCellVolume = 0.;
    Vector<scalar> filterCenter(0.,0.,0.);
    Vector<scalar> currentCellCenter(0., 0., 0.);
    List<Vector<scalar> > filterCenterList;
    List<Vector<scalar> > UBarList;

    volScalarField filter(
        dimensionedScalar("filter", dimensionSet(0,0,0,0,0,0,0), 0.)

    // Build field center list
    for ( int j = 0; j < nY; ++j ) 
        for ( int i = 0; i < nX; ++i)
                    filterXMin + hX * int(i),
                    filterYMin + hY * int(j), 

    // Loop over filter center coordinates
    forAll(filterCenterList, iFilterCenterList)
        // Define filter center
        filterCenter = filterCenterList[iFilterCenterList];
        //Info<< "filterCenter = " << filterCenter << endl;
        // Set filter
        forAll(cellCenterField, iCellCenterField)
            currentCellCenter = cellCenterField[iCellCenterField];
            normInfty = max( mag(currentCellCenter.x() - filterCenter.x()),
                        max( mag(currentCellCenter.y() - filterCenter.y()),
                             mag(currentCellCenter.z() - filterCenter.z()) ) );
            filter[iCellCenterField] = (( normInfty <= filterRadius ) ? 1. : 0.);

        // Compute filter cell volume
        filterCellVolume = fvc::domainIntegrate(filter * cellVolumeField).value();

        // Compute average velocity
        UBarList[iFilterCenterList] = fvc::domainIntegrate(filter * cellVolumeField * U).value() / filterCellVolume;
        //Info<< "UBar = " << UBarList[iFilterCenterList] << endl << endl;
    } // Loop over filter center coordinates

    // Output the UBar field to a custom-format ASCII file 
    // (please note that the call for sqrt is ambiguous due to the damn 
    // using namespace statement)
    // TODO (low priority): modify output directory cleanly (current solution is hack)
    // TODO (low priority): add progression indicator
    fileName outputFile("postProcessing/UBar.csv");
    Info<< "outputFile = " << outputFile << endl;

    OFstream os(outputFile);
    os << "\"x\",\"y\",\"z\",\"UBarX\",\"UBarY\",\"UBarZ\",\"UBarMag\"" << endl;

    forAll(filterCenterList, iFilterCenterList)
        os  << filterCenterList[iFilterCenterList][0] << ","
            << filterCenterList[iFilterCenterList][1] << ","
            << filterCenterList[iFilterCenterList][2] << ","
            << UBarList[iFilterCenterList][0] << ","
            << UBarList[iFilterCenterList][1] << ","
            << UBarList[iFilterCenterList][2] << ","
            << Foam::sqrt(  sqr(UBarList[iFilterCenterList][0])
                          + sqr(UBarList[iFilterCenterList][1])
                          + sqr(UBarList[iFilterCenterList][2])
            << endl;
        if ( (iFilterCenterList+1)%nX == 0 ) os << endl; // Useful for gnuplot pm3d
This is pretty slow (very slow, actually) but it works fairly well. See the PDF attached for a Gnuplotted filtered velocity field associated with the previously mentioned case.

Handling multiprocessor cases is critical, and I found the solution there:
volum integral of mag(U)...?

I believe there's a bottleneck at the filter setting operation, which is not parallelized.
Attached Files
File Type: pdf UBar.pdf (61.8 KB, 21 views)

Last edited by leroyv; February 27, 2014 at 13:54.
leroyv is offline   Reply With Quote

Old   March 6, 2014, 13:56
Question Use of LES filters?
Vincent Leroy
Join Date: Jul 2012
Location: Bordeaux, France
Posts: 38
Rep Power: 5
leroyv is on a distinguished road
Dear foamers,

Does anyone know if an LES filter wouldn't be suitable to do this? I can't exactly understand what the simpleFilter, laplaceFilter, and such, do.
leroyv is offline   Reply With Quote

Old   October 9, 2014, 03:16
New Member
Sebastian Bomberg
Join Date: Aug 2012
Location: Munich, Germany
Posts: 12
Rep Power: 5
sebas is on a distinguished road
OK. I see this thread is rather old but still...

I implemented some explicit filtering/moving average, too. It works in parallel but takes ages. And edge effects are also an issue with this filtering stuff.

So here's what I'm doing now:

In k-space (wave number), the frequency response of a moving average filter (rectangular impulse response) is a sinc function:

<phi> = sin(L/2*k)/(L/2*k) phi

for a generic variable phi, filter length L and brackets <> denoting filtering.

you can express this as an infinite series:

sinc(L/2*k) = sum_n [L/2*(i*k)]^(2*n) /[(2*n+1)!]

for n from zero to infinity.

Note that (i*k)^2 is essentially the Laplacian.

Now let's write this in a recursive manner:
At iteration m
<phi>_m = <phi>_(m-1) + (L/2)^2/[2*m*(2*m+1)] Laplacian[<phi>_(m-1)]

This is basically the heat equation with a time-dependent diffusivity (tau is the numerical time step)

d<phi>/dt = (L/2)^2/tau /[2t/tau*(2t/tau+1)] Laplacian[<phi>]

Actually, from the series expansion we get Explicit Euler time integration which is numerically unstable unless Diffusion number is very small.

Although it's mathematically not that sound, I just use implicit time integration here.

Remember this is "pseudo time". You have to run this until some steady state at every real time step. But Laplace Equation is easy to solve and as the diffusivity decreases with time it should converge rather quickly.


This derivation is for 1D. Physically it makes sense, I think. But I'd be glad if someone could comment on the multidimensional problem...

Last edited by sebas; October 9, 2014 at 09:14.
sebas is offline   Reply With Quote

Old   October 9, 2014, 11:01
Vincent Leroy
Join Date: Jul 2012
Location: Bordeaux, France
Posts: 38
Rep Power: 5
leroyv is on a distinguished road
Hi Sebastian,

This way of doing things seems quite elegant to me. However, there are a few things that remain unclear to me. I'm going to do a more detailed development of the calculus, using proper math notations for more clarity. That will help me explain where the things I don't understand are.

We start with a given field \phi we want to get the average of using what I'll call the m_0 top-hat filter of width L:
m_0(x) = 1/L \text{   if   } -\frac{L}{2} \leq x \leq \frac{L}{2}
and 0 otherwise. Note that this filter is the first of a family of recursively defined filters that have interesting properties regarding control of the regularity of the resulting averaged field (see Angeli et al., 2013).

We get the average field \left< \phi \right> by applying the m_0 filter using a convolution operator. \left< \phi \right>'s expression in x is:
\left< \phi \right>(x) = (m_0 \star \phi)(x) = \int_{-\infty}^{+\infty} m_0(x-t) \phi(t) \mathrm{d}t

Now, we go to the phase space by applying the Fourier transform operator \mathcal{F}:
\mathcal{F}[\left< \phi \right>](k) = \mathcal{F}[m_0](k) \   \mathcal{F}[\phi](k)
where k is the wave vector.

The only problem-independant thing we have here is the Fourier transform of the filter. As you mentioned, the Fourier transform of m_0 is a sinc:
\mathcal{F}[m_0](k) = \mathrm{sinc}\left(\frac{kL}{2}\right)
and its Taylor series expansion is
\mathrm{sinc}\left(\frac{kL}{2}\right) = \sum_{n = 0}^{+\infty} \frac{(ik\frac{L}{2})^{2n}}{(2n+1)!}

This gives us the following expression for \mathcal{F}[\left< \phi \right>]:
\mathcal{F}[\left< \phi \right>](k) = \sum_{n = 0}^{+\infty} \left( \frac{(ik\frac{L}{2})^{2n}}{(2n+1)!} \mathcal{F}[\phi](k) \right) = \lim_{n \rightarrow +\infty} S_n
S_n = \sum_{m=0}^n \frac{(ik\frac{L}{2})^{2n}}{(2n+1)!} \mathcal{F}[\phi](k)

For there, we can move to a recursive definition of the (S_n(k)) sequence:
S_n(k) = S_{n-1}(k) + \frac{(ik\frac{L}{2})^{2n}}{(2n+1)!} \mathcal{F}[\phi](k)
with the initial condition
S_0(k) = \mathcal{F}[\phi](k)

Finding the limit of the (S_n(k)) sequence would clearly give us \mathcal{F}[\left< \phi \right>], indeed. What I don't understand at this point is how you proceed to the laplacian-based formulation you provide: the (ik)^2 product meaning "apply laplacian" only works if you go back to the real space, in which the recursive definition doesn't work. Plus, we don't have (ik)^2 in the recursive definition; only (ik)^{2n}, a lot less convenient. Did I miss something?

I've also been trying to write this for a 2D problem, but I couldn't figure out how to get to a recursive expression for the Taylor series terms, mostly because the Fourier transform of the 2D m_0 filter is a product of two Taylor series (I'm not a good mathematician).

Last edited by leroyv; October 9, 2014 at 11:02. Reason: grammar
leroyv is offline   Reply With Quote


filter, large eddy simulation, les, openfoam, volume averaging

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
No layers in a small gap bobburnquist OpenFOAM Native Meshers: snappyHexMesh and Others 6 August 26, 2015 09:38
draw any scalar vs. cell volume wersoe OpenFOAM Post-Processing 0 April 3, 2012 03:52
mesh airfoil NACA0012 anand_30 OpenFOAM Meshing & Mesh Conversion 12 December 12, 2011 05:16
negative cell volume in dynamic mesh WU zhonghua FLUENT 0 July 28, 2004 10:04
calculate cell volume, center...? Paul Main CFD Forum 5 June 21, 2003 12:55

All times are GMT -4. The time now is 02:46.