Finite difference formulation of the heat equation with thermal conductivity in 1D

Question

This may seem trivial, but I'm having some trouble deriving the finite difference form of the heat equation with a thermal conductivity function $a(x)$ depending on $x$:

$$\frac{\partial u(x, t)}{\partial t} = \frac{\partial }{\partial x}[a(x)\frac{\partial u(x, t)}{\partial x}]$$

I find some trouble deriving the right hand side. Let $h$ be the space step, and $x_{i−1}=x_i−h$ and $x_{i+1}=x_i+h$. I want to evaluate the right hand side on node $i$. If I take centered differences I get:

$$\begin{align*}\left\{\frac{\partial}{\partial x}\left[ a(x)\frac{\partial \, u(x, t)}{\partial x}\right]\right\}_i &= \frac{\left[a(x)\frac{\partial\, u(x, t)}{\partial x}\right]_{i+1/2} - \left[a(x)\frac{\partial \, u(x, t)}{\partial x}\right]_{i-1/2}}{h} \\ &= \frac{a_{i+1/2}\left[\frac{u_{i+1}-u_{i}}{h}\right] - a_{i-1/2}\left[\frac{u_{i}-u_{i-1}}{h}\right]}{h} \end{align*}$$

Is this ok? The conductivities should be in between the nodes I'm calculating? (intuitively it seems this is reasonable). Is there a way, for example by using forward or backward differences, to use conductivity nodes that overlap with the nodes I'm calculating (I calculate only the interger indexed nodes)? Can I mix a centered difference approach for one derivative and a backward or forward approach for the other one?

I would really appreciate any hint.

Thanks in advance,

Federico

Answer 1

At first glance this looks ok. It is also reasonable to calculate your conductivity at cell nodes.

This is a diffusion equation so the only scheme of any real practical use here is central differencing as by the very nature of diffusion the 'flow' of information is isotropic. For heat transfer through fluids (moving fluids) then upwind differencing can come into play to describe the fluid advection, but the difference scheme used for the diffusion terms in these equations would still be central differencing.

For a comprehensive derivation of this scheme and others like it for thermal heat transfer through solid and fluid media, see:

1. Patankar, Suhas V. (1980), Numerical Heat Transfer and Fluid Flow, Hemisphere.

2. Tannehill, John C., et al., (1997), Computational Fluid mechanics and Heat Transfer, 2nd Ed., Taylor and Francis.

3. Versteeg, H. K., Malalasekera, W., (2007), An Introduction to Computational Fluid Dynamics, The Finite Volume Method. Second Edition.

I hope this was of some use.

Answer 2

I would use the chain rule to get $a$ separate from $\partial u/\partial x$, which will end up as $\partial^2/\partial x^2$. To do the second derivative, your method worked fine.

Answer 3

Let us ignore time $t$ and assume that the functions $u:\mathbb{R}\to \mathbb{R}$ and $a:\mathbb{R}\to \mathbb{R}$ are sufficiently smooth functions.

In a discretization of a (higher-order) differential quotient, it is in principle ok to use any (higher-order) difference quotient that reproduces the (higher-order) derivative in the limit of vanishing discretization $h\to0$. OP's proposed discretization

$$ \frac{a(x+\frac{h}{2})\frac{u(x+h)-u(x)}{h} - a(x-\frac{h}{2})\frac{u(x)-u(x-h)}{h}}{h}$$

has the added benefit that it approaches the derivative $$ \frac{d}{d x}\left[ a(x)\frac{d u(x)}{d x}\right] $$ a bit faster for $h\to0$, namely as ${\cal O}(h^2)$.