# How to troubleshoot numerical instability using finite difference for steady-state non-linear heat conduction equation

I have a problem which I believe is numerical instability when trying to solve a heat conduction equation using finite difference. The short version is that when the parameter $$I=80.3$$ I get the blue curve below (positive values and with one local maximum, which is physically reasonable), but when $$I=80.4$$ I get the orange curve which is not physically correct. I would appreciate any help troubleshooting this problem.

Here is the long description, which I divided into physical modelling (for completeness's sake), finite difference discretization and numerical implementation (where I believe lies the problem):

1: Physical Modelling.

The problem I am trying to solve is of a thin straight metallic wire, whose both ends are held at constant temperatures $$T_0$$ and $$T_L$$. Through the wire flows a current $$I$$ which heats up the wire through Joule heating. The electrical resistivity is temperature-dependent as $$\rho(T) = \rho_0 + \rho_1 \cdot T$$

2: Finite Difference Discretization.

I divide the wire in $$Nx$$ nodes, for each inner node I write the following balance equation

$$q_{i-1} + q_{i+1} + q_e = 0$$

Where $$q_{i-1}$$ is the heat coming from node $$i-1$$, $$q_{i+1}$$ is the heat coming from node $$i+1$$ and $$q_e$$ the heat generated by Joule heating.

I apply Fourier law $$q = -k \cdot A \cdot \Delta T / \Delta x$$, where $$k$$ is heat conductivity, $$A$$ is area, $$\Delta T$$ temperature difference and $$\Delta x$$ is the length between two nodes.

$$\frac{k \cdot A}{\Delta x} \cdot (T_{i-1} - T_i) + \frac{k \cdot A}{\Delta x} \cdot (T_{i+1} - T_i) + (\rho_0 + \rho_1 \cdot T_i) \cdot \frac{\Delta x}{A} \cdot I^2 = 0$$

Reordering I get

$$\alpha \cdot T_{i-1} + (-2\alpha + \rho_1 \cdot \frac{\Delta x}{A} \cdot I^2) \cdot T_i + \alpha \cdot T_{i+1} = -\rho_0 \cdot \frac{\Delta x}{A} \cdot I^2$$

3: Numerical Implementation.

I put the left side of the previous equation in a sparse tridiagonal matrix $$A$$, the right side in a vector $$b$$ and use scipy to solve the system of equations using float64. Here is the Python code:

# Imports
import numpy as np
import scipy.sparse.linalg as la
import scipy.sparse as sp
import matplotlib.pyplot as plt

#% Geometric specifications

D_L = 0.4E-3 # in m
L_L = 12E-3  # in m
T_0 = 87  # in °C
T_L = 175 # in °C

# Material constants

rho0 = 1.6E-8 # Ohm m
rho1 = 6.6E-11 # Ohm m / K
k  = 3.94E2 # W/m K

Nx = 30 # Number of intervals

x   = np.linspace(0, L_L, Nx+1)
d_x = x[1] - x[0]

for I in np.array( [80.3, 80.4] ):

u = np.zeros(Nx+1)

# Representation of sparse matrix and right-hand side
main  = np.zeros(Nx+1)
lower = np.zeros(Nx)
upper = np.zeros(Nx)
b     = np.ones(Nx+1)

# Construction of matrix equations

A_c   = 3.14*(D_L)**2 /4
alpha = k * A_c / d_x
R     = rho0 * d_x / A_c
Pe    = R * I**2
b     = -Pe * b

# Precompute sparse matrix
main[:]  = -2*alpha + rho1 * d_x / A_c * I**2
lower[:] = alpha
upper[:] = alpha

# boundary conditions
main[0]  = 1; upper[0] = 0; b[0]  = T_0
main[-1] = 1; lower[-1]= 0; b[-1] = T_L

A = sp.diags( diagonals=[main, lower, upper], offsets=[0, -1, 1],
shape=(Nx+1, Nx+1), format='csc' )

u[:]   = la.spsolve(A, b)

plt.plot( x, u, label='I='+str(I) )

plt.legend()


Does somebody know how to troubleshooting this problem? Many thanks in advance!

• Should I**2 be in main[:] = -2*alpha + rho1 * d_x / A_c * I**2? From the discussion above, it appears that I is only present in b Dec 2, 2021 at 17:24
• Thanks Charlie for your comment, yes it should be there, it was a typo in the Latex code, now it's fixed Dec 2, 2021 at 18:04
• Just a quick comment before lunch ;-) - check the signs in your formulation, simply changing one line in your code from: main[:] = -2*alpha + rho1 * d_x / A_c * I2 to main[:] = -2*alpha - rho1 * d_x / A_c * I2 gives expected results. Dec 3, 2021 at 10:18
• Hi @PeterFrolkovič thanks for you comment, I checked my derivation and the plus sign before rho1 is correct, otherwise it would mean that electric resistivity decreases which temperature. Dec 3, 2021 at 13:03
• Then check the sign before diffusion, typically it results in a positive number on the diagonal main[:], you have the opposite case. Dec 3, 2021 at 13:23