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

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?
• If you have a no-flux boundary, then I think the gradient is zero by definition. If you have an in/out-flux boundary condition, the gradient can be calculated from the prescribed flux. – MPIchael Aug 25 '20 at 9:54
• @MPIchael "If you have a no-flux boundary, then I think the gradient is zero by definition." Yes, you're correct, but consider an unstructured mesh where the face normal vector does not coincide with the vector joining two adjacent cells centroids, the flux - due to non-orthogonality- in this case is decomposed into two parts: 1) Linearized part, where we calculate the gradient in term of $\phi$ at cells centroids divided by the distance between them. 2) Non-linearized part (cross diffusion term), where we compute the gradient through Gauss-Green approximation. – Algo Aug 25 '20 at 10:34
• @MPIchael And my question is mainly about that second part, how to calculate the average gradient of $\phi$ at cell centroid in terms of value of $\phi$ at its faces. In case one of those faces is insulated, what is the value of $\phi$ at that face?. The second part of my question is when I need to calculate the gradient at a boundary face, that is part of only one cell and I cannot perform interpolation. – Algo Aug 25 '20 at 10:37
• In your discretization, what lokal base do you use (Lagrange Base, piecewise constant)? – MPIchael Aug 25 '20 at 14:02
• No-flux boundary condition in general does not mean zero gradient. It depends on the transport model. If transport is diffusive $j = -D \nabla \phi$ that's one thing. But if $j = V \phi$, where $V$ is a given prescribed velocity field, that's something else. Prescribed velocity field can be realized for example for passive advection of particles suspended in a fluid flow. – Maxim Umansky Aug 26 '20 at 4:07

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

• Thanks for following up! – MPIchael Sep 1 '20 at 12:29
• @MPIchael Thanks :)) – Algo Sep 1 '20 at 13:20