Instability in 2D steady state heat equation with variable thermal diffusivity

Question

Thank you in advance for your time and consideration with this issue. Any suggestions would be greatly appreciated!

I'm trying to numerically solve the steady state heat equation in 2D (x,y) with variable thermal diffusivity. I started with the following equation:

$$ 0 = \mathbf{\nabla} \cdot \left( \alpha \mathbf{\nabla}U\right) $$

$$ 0 = \alpha \nabla^{2}U + \nabla\alpha\cdot\nabla U $$

$$ 0 = \alpha \left[ \frac{\partial^{2}U}{\partial x^{2}} + \frac{\partial^{2}U}{\partial y^{2}} \right] + \frac{\partial\alpha}{\partial x} \frac{\partial U}{\partial x} + \frac{\partial\alpha}{\partial y} \frac{\partial U}{\partial y}$$

Where $U$ is the temperature and $\alpha$ is the thermal diffusivity, both are functions of x and y.

I used a centered difference scheme for both the temperature $U$, and the thermal diffusivity $\alpha$. Solving for $U_{i,j}$ I get:

$$ U_{i,j} = \frac{1}{4} \left[\left(1 + \frac{\alpha_{x}}{2\alpha} \right)U_{i+1,j}+\left(1-\frac{\alpha_{x}}{2\alpha}\right)U_{i-1,j} +\left(1+\frac{\alpha_{y}}{2\alpha}\right)U_{i,j+1} + \left(1 - \frac{\alpha_{y}}{2\alpha}\right)U_{i,j-1} \right]$$

Where the finite difference forms of the first derivative of thermal diffusivity are:

$$ \alpha_{x} = \frac{\alpha_{i+1, j} - \alpha_{i-1, j}}{2} $$

$$ \alpha_{y} = \frac{\alpha_{i, j+1} - \alpha_{i, j-1}}{2} $$

I initialize a 2D grid with a starting temperature (100) and set a boundary condition on one side (200). As I iteratively solve for $U$ many of the initial values of $U$ drop below 100 (99.99999999999997) during the first iteration and eventually the iterations diverge (i.e. become negative).

Before I show the code that I wrote to implement this here is a list of things I've tried:

• Modified the thermal diffusivity matrix to have all the same values (i.e. constant thermal diffusivity). This actually worked and the solution converged. As one expects the derivative terms of thermal diffusivity are zero. This makes me think that it has something to do with the finite difference representation of the thermal diffusivity.
• Created all the variables using double precision (128). The iterations still became unstable.
• I tried both rescaling the thermal diffusivity by a single value and also a linear transformation to decrease the range of values. The iterations still became unstable.

Here is the snippet of the Python code I wrote to solve this:

import numpy as np

def partial_k(k):
    k_x = 0.5 * (k[2:, 1:-1] - k[:-2, 1:-1])
    k_y = 0.5 * (k[1:-1, 2:] - k[1:-1, :-2])

    return k_x, k_y


def delta(u, k, k_x, k_y):
    # Propagate with central-difference in space
    u_new = u.copy()
    u_new[1:-1, 1:-1] = 0.25 * (
        (1 + 0.5 * k_x / k[1:-1, 1:-1]) * u[2:, 1:-1] +
        (1 - 0.5 * k_x / k[1:-1, 1:-1]) * u[:-2, 1:-1] +
        (1 + 0.5 * k_y / k[1:-1, 1:-1]) * u[1:-1, 2:] +
        (1 - 0.5 * k_y / k[1:-1, 1:-1]) * u[1:-1, :-2]
    )

    return u_new


# Files and directories
fname = 'some/file/with/thermal/values'

iter_out = 500

# Get Friction surface from raster
k = np.load(fname)

# Find partials of thermal conductivity
k_x, k_y = partial_k(k)

# Number of rows and columns
n_rows = np.shape(k)[0]
n_cols = np.shape(k)[1]

# Max number of iterations
max_iters = 5001

# Convergence tolerance
tol = 1e-5

# Set the boundary conditions
t_cool = 100
t_hot = 200

# Initialize temperature matrix
u_0 = t_cool * np.ones((n_rows, n_cols))

# Set boundary conditions
u_0[:, -1] = t_hot

# -----------
# Calculation
# -----------
for i in range(max_iters):

    u_update = delta(u_0, k, k_x, k_y)

    delta_u = np.abs(u_0 - u_update)

    if delta_u.max() < tol:
        break
    else:
        u_0 = u_update.copy()