Flux sign and face normal confusion in finite volume method

An extended answer.

For more arbitrary meshes you have to consider that generally CFD/FEM solvers rely on generic data-structures with element and side lists:

Element list

Side list

Consider the following pictures, which is the standard case for simple Cartesian meshes. Since there is a single plus and a single minus side on each face, the definition is unique and calculations does not cause any problems. In your words this would mean

• East sides always match west sides.
• South sides always match north sides.

However, due to performance reasons and due to parallelization strategies fluxes and derivatives at the faces are generally not calculated element wise. A naive implementation would result in unnecessary double calculations. A simple definition via plus and minus no longer unique.

For generic data structures you have to consider two things:

1. Elements may be rotated between each other.
2. Elements may have a different orientation.

This has the consequence that two same sides would match each other. Due to this, it is necessary to define an additional mapping to get a unique side list, here defined as master and slave sides. Moreover, due to a wrong orientation tangential vectors may be pointing in a different direction. This may become important for high order elements. The following pictures contain elements rotated between each other. Here each master side matches a slave side.

When dealing with conservation laws like your case, you can often make use of the divergence theorem (as you did). You can then express the fact that the total mass within your integration region is preserved by the following surface integral:

$$oint_{partial Omega} k nabla T cdot mathbf{n} ~partial S = 0$$

Now, as it stands, it is irrelevant which direction $mathbf{n}$ is pointing to, as long as it is consistent across this surface integral.

When you take the convention that the direction of $mathbf{n}$ is outward pointing, you say that the total mass being transported out is zero. When you say the normals point to the inward, you integrate the flux into the cell, so to say. The statements are equivalent. Proper care has to be taken to have the signs consistent when you have source terms.

In more mathematical language, you may have two perfectly correct formulations of your gauss theorem, one with inward pointing normals, and one with outward pointing normals:

canonical, outward pointing gauss-theorem: $$ int_V operatorname{div} vec F ; mathrm d^{(n)}V = oint_{S} vec F cdot vec n_{out}; mathrm d^{(n-1)}S$$

noncanonical, inward pointing gauss-theorem: $$ int_V operatorname{div} vec F ; mathrm d^{(n)}V = - oint_{S} vec F cdot vec n_{in}; mathrm d^{(n-1)}S$$

The statements are of course equivalent. In the second case you just have a sign on the r.h.s.

I do not want to add to the confusion. The consensus is, to always have outward pointing normals for any given control volume, and many conservation laws are written in that way, but there is an element of convention to it. If you employ this convention in your code consistently, the mass flowing outside of one cell (negative sign) will be added to the adjacent cell.

Normal direction depends on the cell that you are writing equation for. the word outward is relative to the cell under study. In order to write equation for each of cells, i.e. $Sigma nabla T.n S_f=0$, stick to this : $nabla T_{face}=frac{T_c-T_i}{r_c-r_i}$ and assume $n$ as outward pointing normal vector for that face. I think your problem is that you have mistakenly introduced two negative signs in some places of your code; one in the denominator of $frac{T_c-T_i}{r_c-r_i}$ and one in the normal vector $n$. OpenFOAM classifies the cells as owner and neighbour to determine outward direction.