Flux sign and face normal confusion in finite volume method

Question

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)

1. 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?

2. 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.

Answer 1

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.

Answer 2

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.

Answer 3

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.