How to solve the heat equation for compound materials with different heat conductivities numerically?

Question

I'm solving the heat equation with time dependent boundary conditions numerically in a 2D system using the ADI scheme. For the purpose of this question, let's assume a constant heat conductivity and assume a 1D system, so $$ \rho c_p \frac{\partial T}{\partial t} = \lambda \frac{\partial^2 T}{\partial x^2}. $$ This works very well, but now I'm trying to introduce a second material. This one differs slightly in heat capacity and density but has a very different heat conductivity and is connected to the other material by a sharp interface, i.e. a stepwise change in $\lambda$.

How should this be treated numerically in the ADI scheme? I can think of different approaches:

1. Treat the two materials as independent domains and connect them by a boundary condition calculating the heat flow in and out of the interface in terms of temperature on the other side of the interface in the last time step. Use a simple forward difference for that on both sides of the interface.
2. Treat it as one domain and use a very fine discretization close the interface as compared to the homogeneous material. Use a scheme like $$ \lambda_{left} \frac{T_i - T_{i-1}}{\Delta x_i} = \lambda_{right} \frac{T_j - T_{j+1}}{\Delta x_j}, $$ where $i$ and $j$ are the points left and right of the interface, instead of the the standard ADI for those points.
3. Drop the assumption of constant heat conductivity and use $$ \rho c_p \frac{\partial T}{\partial t} = \frac{\partial \lambda}{\partial x}\frac{\partial T}{\partial x} + \lambda \frac{\partial^2 T}{\partial x^2}. $$ But in order do so one needs to approximate derivative of lambda at the step position, i.e. introduce an unknown characteristic width $s$ of the sharp interface. I assume, that the (more or less arbitrary) choice of this width will significantly influence the system's behaviour.

Any advice?

Answer 1

The generalized heat equation is, $$ \frac{\partial T}{\partial t}=\nabla\cdot\left(\alpha\nabla T\right) $$ If $\alpha$ is spatially independent, then we can pull it out of the differential operator and obtain your 1st equation. If $\alpha$ is spatially dependent, then numerically in one dimension, you have $$ \frac{T^{n+1}_i-T^n_i}{dt}=\frac1{dx^2}\left[\alpha_{i+\frac12}\left(T^n_{i+1}-T^n_{i}\right)-\alpha_{i-\frac12}\left(T^n_i-T^n_{i-1}\right)\right] $$ where $$ \alpha_{i+\frac12}=\frac12\left(\alpha_{i+1}+\alpha_{i}\right) $$ and all other terms take their normal meaning. Extension to 2D or implicit methods should be trivial from here.