Computing face fluxes in FVM

Question

In FVM, we have to compute fluxes at some face of a cell. There are many ways to compute this face flux value, but the most common and easiest way involves some simple averaging at the face. However, I am a little confused on what the best way of averaging is.

The particular case I am wondering about is when you have some material property multiplied by a gradient. For example, the heat conduction equation is $k \frac{dT}{dx}$.

For this, we want to compute $k \frac{dT}{dx}$ at a face. By using the simple averaging technique, I can see two ways you can do this.

To illustrate the two ways, let's assume that the cell $L$ and cell $R$ are the neighbors of the face. Let's assume the cell centers are of equal distance from the face centroid, to avoid talking about weighted averages.

The two ways are:

(1) $$ 0.5(k_L\frac{dT}{dx}|_L + k_R\frac{dT}{dx}|_R) $$

(2) $$ 0.5(k_L + k_R) * 0.5(\frac{dT}{dx}|_L + \frac{dT}{dx}|_R) $$

In the first case, we compute the heat conduction in each cell and then average it to get the heat conduction at the face. In the second case, we average the conductivity in each cell, and also the temperature gradient in each cell and then we compute the heat conduction at the face from the averaged conductivity and average temperature gradient.

Is one method superior vs. the other? Obviously, both methods are equivalent if $k_L = k_R$, so this discussion assumes they are different.

Answer 1

You are looking for a harmonic average. If you think about the face between two cells, what kind of peculiar cases can happen?

1. If your material properties are equal, you want the flux to represent that.

2. If one of your materials is impermeable (very low heat conduction in your case) you want your flux to be zero!

3. In any other case you want to use a reasonable value between the material properties in the adjacent cells.

The canonical way is, I think, to use the harmonic average:

$$k_\text{interface} = \frac{1}{(1/k_l) + (1/k_r)}$$

If you squint at that averaging, you see that a $k_l = 0$ or a $k_r=0$ will force your interface flux to be zero. If the values are well above zero, but different, then the average will give a reasonable value between the two. If they are equal, you will also have the right average.

It gets a little more advanced if your material properties are not scalar, i.e., vectors or, god forbid, diffusion tensors, but in your case that should suffice.

For the derivatives you need you have to consider how many degrees of freedom you have in your cell. If you have a structured quadrilateral grid with cell centered DOFs, then you would discretize the derivative over the interface with the two dofs in the adjacent cells. If you happen to have more dof's per cell you have more choices, but the canonical FV way is with one dof per cell.

So in total you would have something like this:

$f = \frac{u_r - u_l}{h} ~ k_\text{interface} ~A$

With the discretized derivative, the averaged material property $k_\text{interface}$, and the face area A.

If you want to dig in deeper, on why that harmonic average is the way to go check out the intro in chapter 3 ("Interaction volume") in the following link. You can propably skip most of the paper, but they may put it better than I can: Ivar Aavatsmark - Multipoint flux approximation methods forquadrilateral grids