Support for two term FV gradient expansion with triangle mesh

[GiudGiud]

Bug Description

In INSFV simulations:
When using triangle mesh with a mesh cell in a corner, which has both a wall and a flow boundary condition, the matrix (the one solved to obtain the face gradients, the centroid gradient and the face values of velocity) is reported as singular and the simulation diverges.
In general FV simulations, it's likely the issue will appear with future ZeroNormalGradient boundary conditions

Steps to Reproduce

Any mesh with a triangle cell that covers an entire corner. It's likely though unconfirmed that splitting the corner with two mesh cells is sufficient to remove the problem. Indeed if the corner is split as is done automatically by gmsh (see below) then there is no issue. The diverging channel tests in INSFV use the below mesh with two term boundary expansion.

It's likely though unconfirmed the issue can also be obtained on quad meshes with different combinations of boundary conditions (3 boundaries for a cell for example). It has been confirmed that if three faces of a (orthogonal) quad are on a boundary, then singularity occurs

Impact

This forces the use of quad meshes in INSFV simulations.

[lindsayad]

Yea I don't think there's much I can do about this. I worked this out by hand for the tri test case in the postprocessors directory, and the matrix is just straight up singular. I could resolve it by projecting dCb onto the surface normal vector ... that would guarantee to remove the singularity. Essentially then you are just extrapolating the boundary face value with the part of the gradient that is in the face normal direction

[lindsayad]

This patch allows me to run the postprocessor tri test case with two term boundary expansion. If you think this is a satisfactory enough solution, then I will submit a PR

diff --git a/framework/src/variables/MooseVariableFV.C b/framework/src/variables/MooseVariableFV.C
index 11bee5484c..4c29825a62 100644
--- a/framework/src/variables/MooseVariableFV.C
+++ b/framework/src/variables/MooseVariableFV.C
@@ -782,9 +782,15 @@ MooseVariableFV<OutputType>::adGradSln(const Elem * const elem) const
         ebf_faces.push_back(fi);
 
         // eqn. 2
-        ebf_grad_coeffs.push_back(-1. * (elem_has_info
-                                             ? (fi->faceCentroid() - fi->elemCentroid())
-                                             : (fi->faceCentroid() - fi->neighborCentroid())));
+        const auto surface_normal = surface_vector / surface_vector.norm();
+        const auto dCb = elem_has_info ? (fi->faceCentroid() - fi->elemCentroid())
+                                       : (fi->faceCentroid() - fi->neighborCentroid());
+        // We project dCb onto the surface_normal direction in order to prevent singularities when
+        // we have multiple extrapolated boundary faces per element
+        mooseAssert(dCb * surface_normal >= 0,
+                    "I think these things should be in the same direction...");
+        const auto dCbn = (dCb * surface_normal) * surface_normal;
+        ebf_grad_coeffs.push_back(-1. * dCbn);
         ebf_b.push_back(&elem_value);
 
         // eqn. 1

[GiudGiud]

I like that, but I want to know more

Does that break all the other tests?

Do you think it's the same scenario for the extrapolation with tetrahedra?

Is it only when two of the three faces are extrapolated? Or already with a single extrapolated face?

I can't see anything in Moukalled about a special treatment for those meshes. Can you think of another reference that might have that?

[lindsayad]

I like that, but I want to know more

🤣

Does that break all the other tests?

About to find out

Do you think it's the same scenario for the extrapolation with tetrahedra?

Is it only when two of the three faces are extrapolated? Or already with a single extrapolated face?

It requires multiple extrapolated faces and it requires that for at least one of the extrapolated faces that (dCf * n) n != dCf

I can't see anything in Moukalled about a special treatment for those meshes. Can you think of another reference that might have that?

Moukalled doesn't seem to go into detail on a lot of the finer points. For instance skew corrections get half a page and there is no mention of solving for a gradient with Green-Gauss when you have extrapolated boundary faces. This article goes into quite a bit of detail on gradient computation but seems to imply that Green-Gauss is second order at internal faces and first order at boundary faces ... which is only true if you're assuming you're going to use the internal cell value as the boundary face value.

[lindsayad]

Or perhaps they're referring to the accuracy of the gradient itself? In that case I'd be shocked if Green-Gauss yields a second-order gradient. Usually if an error norm for the solution is order N, then a similar error norm for the gradient is order N - 1.

[lindsayad]

I ended up deciding to go the route of just using one-term boundary expansion by default for INSFV. You can see my explanation in #17716. If we keep the current title of this issue, then I don't think #17716 closes it, but if we amend the title to "cannot run INSFV with many triangle/tet meshes" then #17716 can be considered a fix. Up to you

[GiudGiud]

I think it's the right fix for us but if we keep it open it may avoid some confusion later on
Mesh convergence studies are pretty much standard procedure for CFD, so I can see the question come up soon enough