I am modifying OpenFVM to deal with the neutron diffusion equation. The equation consists of a second order elliptic operator, in a way identical to the viscous diffusion term in the momentum equation. I noticed that when I perform the non-orthogonal correction, the code "blows up" if I divide the cell-center gradient by the volume of the mesh (Vp). If I don't divide by it, everything seems to be fine, meaning I get convergence and the results "make sense" compared to a case where I do not apply any non-orthogonal correction.
Non-orthogonal corrections tend to make the solution more unstable because it adds "weight" to the source term. With non-orthogonal corrections it might be necessary to reduce the time step significantly.