# How to correctly define the flux in a finite volume method for Poisson's equation with a piecewise constant material

How do we correctly define the flux in a finite volume method applied to Poisson's equation where we have a piecewise constant material? Specifically, say we have the equation \begin{align*} -\nabla\cdot(k\nabla u) =& f\textrm{ on }\Omega\\ u =& 0\textrm{ on }\partial\Omega \end{align*} where $\Omega$ is some kind of easy, square domain in $\mathbb{R}^2$. Then, let's mesh $\Omega$ using square elements where one of the cells is $K$. Then, on the interior, we can integrate over $K$, which gives $$\int_{K}-\nabla\cdot(k\nabla u) = \int_{K}f$$ After using the divergence theorem, we have $$\int_{\partial K}-k(\nabla u\cdot n) = \int_{K}f$$ Now, we have to discretize the quantity on the left, which can be divided into four separate integrals. Assume that center of cell $K$ is $(x_i,y_j)$ and let $u_{ij}=u(x_i,y_j)$ and $k_{ij}=k(x_i,y_j)$. Then, we can discretize $\nabla u \cdot n$ into four different pieces, depending on what side of the square we're on, as \begin{align*} &\frac{u_{ij}-u_{{i-1},j}}{\Delta x}\\ &\frac{u_{ij}-u_{{i+1},j}}{\Delta x}\\ &\frac{u_{ij}-u_{i,{j-1}}}{\Delta y}\\ &\frac{u_{ij}-u_{i,{j+1}}}{\Delta y} \end{align*} What I'm not clear about is what to do about $k$ if we assume that $k$ is piecewise constant and defined on the same mesh as $u$. We could just multiply the fluxes by $k_{ij}$. For example, discretize $-k(\nabla u\cdot n)$ along the left face as $$k_{ij}\left(\frac{u_{ij}-u_{{i-1},j}}{\Delta x}\right)$$ However, in this case, the flux on the other side would be $$k_{{i-1},j}\left(\frac{u_{{i-1},j}-u_{i,j}}{\Delta x}\right)$$ so they wouldn't match and be conservative because they're using two different material values. The best I can think of is to average the material and use $$\left(\frac{k_{ij}+k_{{i-1},j}}{2}\right)\left(\frac{u_{ij}-u_{{i-1},j}}{\Delta x}\right)$$ but I'm not sure if this is correct or the implications of this choice.

As such, how do we discretize the flux when we have a piecewise constant material, $k$, defined on the same mesh as our material $u$.

## 1 Answer

Taking the average will certainly work, but is recommended to take the harmonic mean to account for the changing conductivity:

$k = \Big( \frac{1-f}{k_i} + \frac{f}{k_{i+1}} \Big)^{-1}$

If you choose your faces directly midway between your two cells, f = 0.5:

$k = \frac{2k_ik_{i+1}}{k_i + k_{i+1}}$

It can be derived by looking at the physical flux, which you want to approximate.

I suggest, you take a look at the book by Patankar: Numercial Heat Transfer and Fluid Flow Chapter 3. It can be found in any decent science library and I guess, the internet counts as well.