How to compute gradient of a cell having a boundary face?

Question

In many situations in unstructured mesh solvers, one needs to compute gradient of arbitrary variable $\phi$ such as temperature or velocity at face centers (one of such situations is correction for mesh non-orthogonality), which can be approximated as:

$$ \nabla \phi_f = (g_C \ \nabla \phi_C) + (g_F \ \nabla \phi_F)$$

Where $C$ & $F$ are two adjacent cells sharing the face $f$ and $g_C$ & $g_F$ are interpolation weight factors.

Now, $\nabla \phi_C$ and $\nabla \phi_F$ can be approximated as the average gradient across the cell, using Gauss-Green theorem as: $$ \nabla \phi_C = \frac{1}{V_C} \sum_{faces} \phi_f \mathbf{S_f} $$

where $\mathbf{S_f}$ is the face normal vector at face $f$, and $V_C$ is the cell volume.

The $\phi_f$ values for interior faces can be easily known.

1. What if the face is a boundary face and insulated? how to get $\phi_f$ in such case?
2. What if the face has a fixed boundary condition $\phi_{f} = \phi_{specified}$ (meaning that face has only one adjacent cell, so first interpolation equation cannot be used), how to compute the gradient at the face in such case?

Answer 1

Disclaimer: I'm not 100% sure but I thought that I should provide my working solution to the above problem, for any future visitor that might have the same questions. The answer is for a steady temperature diffusion problem.

What if the face is a boundary face and insulated? how to get $\phi_f$ in such case?

Since a boundary face is in contact with only one cell, we cannot use interpolation, but since an insulated boundary will have a zero gradient $\nabla \phi=0$ across all its faces, and by assuming linear variation between cell centroid and face centroid, thus for an arbitrary cell with a face boundary one can write: $$ \phi_C = \phi_f + (r_C - r_f).\nabla \phi_f \implies \phi_C = \phi_f$$

What if the face has a fixed boundary condition $\phi_{f} = \phi_{specified}$ (meaning that face has only one adjacent cell, so first interpolation equation cannot be used), how to compute thegradient at the face in such case?

Again, since a linear profile is assumed: $$ \nabla \phi_f = \frac{\phi_{specified} - \phi_C}{\|r_f - r_C\|} \hat{e}$$ where $\hat{e} = \frac{r_f - r_C}{\|r_f - r_C\|}$

I tried this on an unstructured non-orthogonal mesh and the solution finally converged, having a good agreement with laplacianFoam (I still didn't implement correction for skewness, which justifies the deviation from OpenFOAM (and not enough iterations), and used simple Gauss-Seidel solver), and this is the center line temperature distribution

Answer 2

I had the same question a few days ago when I was developing an unstructured FVM solver for conduction/diffusion. And now I find that no matter what the boundary condition is, just set the gradient of the boundary face the same value as its owner cell's:

$$ \nabla \phi_f = \nabla \phi_C $$

Even for insulated boundary, the gradient of boundary face can not be set to zero, because zero gradient means the cross diffusion can not be corrected for non-orthogonal grid.

My FVM solver is here: https://github.com/xu-xianghua/knablat, you can check the boundary gradient in function correct_cross_diffusion in solver.py .