I implemented a solver for the 2D steady-state heat equation (without heat generation and homogeneous material) $\nabla. (k\nabla T) = 0$, using finite volume method, however, I am having some confusion about flux direction and face normal.
First, using divergence theorem (applied to a mesh element): $$\int_V \nabla. (k\nabla T) \ dV = \oint_S (k\nabla T .n) \ dS = 0$$
Where $n$ is the outward pointing unit vector normal to $S$. And the equation can be further discretized to (using Gaussian quadrature with single integration point):
$$ \sum_{faces} \ (\nabla T .n) S_f = 0$$
My Python solver imports a simple Cartesian (100 x 100 x 1) OpenFOAM mesh (boundary, points, faces, owner and neighbour), applies the discretized equation to each cell and generate the sparse coeffecient matrix $A$ such that $AT = b$.
Originally, I had wrong results because I found out that I was doing the following:
Each interior cell $c$ has four adjacent cells ($n$: north, $s$: south, $e$: east, $w$: west)
- When expressing $\frac{\partial T}{\partial x}$ by assuming that $T$ varies linearly between cells centroids ($c$ and arbitrary adjacent cell $i$): $$\frac{\partial T}{\partial x} = \frac{T_c - T_i}{d_{ci}}$$
First error: I used this formula for all the adjacent cells, which is correct for east and north cells, but when evaluating west and south faces it should be $T_i - T_c$.
First question: I thought that the order of the subtraction won't matter since it will be corrected by the flux direction (next question). So why does the order of the subtraction matters?
- After importing the owner-neighbor OpenFOAM mesh relations, I assumed that the flux is always in the direction from owner to neighbor, so $n$ basically was pointing out of the cell for the faces owned by the cells, and pointing inward for neighbor faces.
Second error: This also resulted in wrong coefficient matrix, and when I corrected $n$ to be always pointing out of the cell faces everything was okay.
Second question: Should I always make the face normals pointing out of the element? If, so how is this physically correct to assume that fluxes are always leaving the cells? (And why having owner-neighbor relationships in the first place).
Sorry, for the long question but I think by providing my full methodology, my confusion will be clear to the reader.