Discretize derivative of heat flux with variable conductivity

Question

How do i discretize this Temperature flow equation with the Fnite Difference method, when the conductivity K, is not constant?: $$ \frac{\partial}{\partial x}\left(K\frac{\partial T}{\partial x}\right) + \frac{\partial}{\partial y}\left(K\frac{\partial T}{\partial y}\right) = 0 $$

My attempt at this problem was to do the following:

$$ K(x)\frac{\partial^2 T}{\partial x^2} + K(y)\frac{\partial^2 T}{\partial y^2} = 0 $$ then discretize the second order partial equations: $$ K(x) \left(\frac{T_{i+1,j} -2T_{i,j} + T_{i-1,j} }{\Delta x}\right) + K(y) \left(\frac{T_{i,j+1} -2T_{i,j} + T_{i,j-1} }{\Delta y}\right) = 0 $$

However, i have been informed that this is an incorrect method as this working assumes a constant conductivity, where that is not the case in the problem.

would someone be able to help me with this?

Answer 1

Let us consider the 1D case for simplicity. The second-order centered finite-difference may be viewed as the composition of two first-order centered finite-differences: $$ \left.\frac{\text d^2 T}{\text d x^2} \right|_{x=x_i}\simeq \frac{T_{i+1}-2T_{i}+T_{i-1}}{{\Delta x}^2} = \frac{1}{\Delta x}\left(\frac{T_{i+1}-T_{i}}{{\Delta x}} - \frac{T_{i}-T_{i-1}}{{\Delta x}}\right) . $$ In the spatially-varying conductivity case, the latter finite-difference writes as $$ \left.\frac{\text d}{\text d x}\left(K\frac{\text d T}{\text d x}\right) \right|_{x=x_i}\simeq K_{i+1/2}\frac{T_{i+1}-T_{i}}{{\Delta x}^2} - K_{i-1/2}\frac{T_{i}-T_{i-1}}{{\Delta x}^2}\, , $$ where $K_{i\pm 1/2} = K\big(x_i\pm\frac{1}{2}\Delta x\big)$. It is incorrect to take $K$ out of the derivative as proposed in OP, since by the product rule, $(K T')' = K T'' + K'T'$ where prime $'$ denotes differentiation w.r.t. $x$. This expression suggests to compute the finite-difference \begin{aligned} \left.\frac{\text d}{\text d x}\left(K\frac{\text d T}{\text d x}\right) \right|_{x=x_i} &\simeq \frac{K_{i+1/2}+K_{i-1/2}}{2}\frac{T_{i+1}-2T_{i}+T_{i-1}}{{\Delta x}^2}\\ &\phantom{ \simeq } + \frac{K_{i+1/2}-K_{i-1/2}}{{\Delta x}}\frac{T_{i+1}-T_{i-1}}{2\,{\Delta x}}\, \end{aligned} which leads to the same method as before. Other suggestions and the same one can be found here. In 2D, we have $\partial_x(K \partial_x T) + \partial_y(K \partial_y T)$ with conductivity $K(x,y)$ instead, but the same principle applies (see e.g. this post).