Hi all -

I wanted to check if this approach is correct for the heat equation with multiple materials (different values for k, rho, cp). For simplicity, I will assume 1D. I am using an implicit scheme.

I assume I have a uniform grid in space and time, with h=\Delta x.

I assume the interface between the two materials lies directly between two nodes. This means that for the last node in the left material, I create a ghost node that coincides with the first node in the right material; these temperatures are donated TL and T{L h}. Vice verse for the right material: a ghost node coincides with the last node in the left material, with nodal temperatures TR and T{R-h}.

To solve for the ghost node in terms of known values, I do the following:

Create a new node at the interface. Both the left material and the right material should yield the same temperature and the same flux at this node. These constraints will allow us to solve for the ghost nodes.

We can create two equations, one to constrain the temperature at the interface and one to constrain the flux.

The temperature at the interface from the “perspective” of the left material is the average of the last node in the left material and the ghost node. Similarly, the temperature at the interface from the perspective of the right material is the average of the ghost node and the first node in the right material.

T{L h/2} = (T_L T{L h})/2

T{R-h/2} = (T_R T{R-h})/2

The flux at the interface from the perspective of the left material can be computed using a central finite difference. Similarly for the flux at the right material.

kL * delT/delx @ L h/2 = k_L * (T{L h} - T_L) / h

kR * delT/delx @ R-h/2 = k_R * (T_R - T{R-h} / h

We end up with the following two equations:

(TL T{L h})/2 = (TR T{R-h})/2

kL * (T{L h} - TL) / h = k_R * (T_R - T{R-h} / h

We can use these two equations to solve for T{L h} and T{R-h} in terms of k_L, k_R, T_L, T_h, and h.

Once solved for, we can then easily plug into the usual discretization of the heat equation at nodes L and R and determine the coefficients for the tridiagonal matrix.

How does this look? Am I missing something obvious? Are there any accuracy issues I should be aware of?

Follow up question: if I instead wanted to assume “rough” contact with a heat transfer coefficient H, is the approach to say that kL * (T{L h} - T_L)/h = H * (T_R-T_L), solve for the ghost node, and then plug that in (and do the same for the right hand side)? I believe that in this approach, the L and R nodes are actually effectively coincident, but I’m not sure.