[LuComet]

Dear foamers,

I have a 2D turbulent flat plate boundary layer simulated using simpleFoam. The case has x = [-20,22], y = [0,20], z = [0,1], standard boundary conditions (velocity inlet, pressure outlet, symmetry, slip walls and no-slip wall), with flow along x (see attached grid.png).

I am running areaNormalIntegrate on arbitrary planes normal to x, bounded or not, using a post processing function:

x1 5;
x1a #eval {$x1-0.1};
x1b #eval {$x1+0.1};

y0 0;
y1 2;

z0 0;
z1 1;

sampling_post
{
    // Mandatory entries (unmodifiable)
    type            surfaceFieldValue;
    libs            (fieldFunctionObjects);

    // Mandatory entries (runtime modifiable)
    fields          ( U );
    operation       areaNormalIntegrate;

    // Calculate on user-defined bounded plane
    regionType      sampledSurface;
    name            plane_bound;
    sampledSurfaceDict
    {
        type            plane;
        source          cells;
        store           true;
        interpolate     false;
        triangulate     false;
        sampleOnExecute true;

        planeType       pointAndNormal;       
        pointAndNormalDict
        {
            // normal x
            normal (1 0 0);
            point  ($x1 0 0.5);

        }
        // clip y
        //bounds ($x1a $y0 $z0) ($x1b $y1 $z1);

    }

    // Optional entries (runtime modifiable)
    surfaceFormat   vtk;
    writeArea       false;

    // Optional (inherited) entries
    writePrecision  8;
    enabled         true;
    log             true;
    writeFields     true;
    writeToFile     true;

}

Launched on my reconstructed solution:

postProcess -func sampling_post -latestTime

I have two problems:

1) If my plane does not fall on a row of nodes, everything seems ok. Integrating the field U, and considering planes for instance at x = 4.99 and 5.01 (node at x = 5):

x = 4.99 (not on node)

surfaceFieldValue sampling_post write:
    total faces   = 144
    total area    = 20

    areaNormalIntegrate(plane_bound) of U = 115.59995

x = 5.01 (not on node)

surfaceFieldValue sampling_post write:
    total faces   = 144
    total area    = 20

    areaNormalIntegrate(plane_bound) of U = 115.60005

The total area is correctly calculated (the yz section of the case is 20x1) as well as the integral (5.78 m/s * 20 m2 = 115.6 m3/s). However, if the plane falls on a row of nodes, I get unexpected results:

x = 5 (on a node)

surfaceFieldValue sampling_post write:
    total faces   = 159
    total area    = 20.310588

    areaNormalIntegrate(plane_bound) of U = 117.39569

The function detects a different number of faces, hence the integration is affected.
You can see from the attached images not_on_node.png and on_node.png that also the cell normals are modified.
The behaviour is the same for any other row of nodes.

2) I also want to perform this integration on a plane clipped at different y, using the bounds instrucion within the sampledSurfaceDict reported above.

I noticed that the cut ends always on the node closest to the specified bounds. In the two images clip1.png and clip2.png, the blue lines created in ParaView show the expected clipping height y = 2 and y = 2.2 specified through the bounds instruction, which differ from the actual extension of the clipped plane.

I tried OpenFOAM-2012 and OpenFOAM-2112 with no differences. I don't know if it matters, but I created the grid through these steps:
- blockMesh
- snappyHexMesh (without layers)
- extrudeMesh (to get a 2D domain)
- snappyHexMesh (only layers)

Anyone knows an explaination or a workaroud?
Thank you so much!

[olesen]

This is an interesting problem, it might be related to this:https://develop.openfoam.com/Develop.../-/issues/2210 and https://develop.openfoam.com/Develop...e_requests/488

I thought that at least with OpenFOAM-v2112 this would be working correctly, but perhaps there is still something not quite right when the cuts align exactly with grid points. Or it could be that they almost align exactly, but there are still slivers on each side (not really sure). Perhaps for these particular cases you also need a more aggressive (full) face regularisation (see the bottom pics on the merge-request).
Presumably your tet cuts will be clean enough not to generate holes there.

[olesen]

BTW: I just noticed that you have "plane" instead of "cuttingPlane" for the sampling - which means that what I just wrote probably doesn't apply since you are using a cell cutter instead of a iso-surface cutter for the planes. Yes, I know that the naming here is completely horrible and non-intuitive, but the sampled planes started off their life by cutting cells and then "walking" the edge cuts to piece together the result. Generally works well for regular cells, but falls apart quickly for complex polyhedral cells. This was partly reworked with some heuristics to handle the edge cases you have (cuts aligned with cell/cell boundaries), but still remains not particularly amazing.
Later the sampled planes were evolved to "cut" cells as iso-surfaces. This newer version was termed "cuttingPlane" (non-obvious naming). This too has evolved from the original node-based algorithm, to a cell-based algorithm and now we've landed at a topology-based algorithm, which generally performs reasonably well.