I am working on a mixed material model for a melting material and need to enforce both a Dirichlet and Neumann type condition at the interface. Subject to an external surface heat flux at the top of the liquid layer, as the material melts, it starts to form a layer of liquid on top of the solid material. I am trying to model the temperature distributions in each domain.
Currently, I am treating them as two separate domains and trying to find a way to couple the two such that the resulting temperature distributions make physical sense. At the boundary, both domains should match in temperature (both at the melting point of the material) and heat flux. For the Crank-Nicolson approach, I have the tridiagonal matrices for each domain:
$$\begin{bmatrix}B&C&&\cdots&0 \\ A & B & C & &\\&A&B&C&\\ \vdots&&A&B&C\\0&&&A&B\end{bmatrix}\begin{bmatrix}T_0^{n+1}\\\vdots\\T_I^{n+1}\end{bmatrix}=\begin{bmatrix}f_0\\ \vdots\\f_I\end{bmatrix}$$
For the upper domain, I am applying a constant heat flux (setting $B=1, C=-1$ and $f_0=q^{''}\Delta x / \kappa$) and a constant temperature boundary condition on the bottom of the liquid domain ($f_I = T_{melt}$). My issue is how to correctly set the boundary conditions for the top of the solid domain. The heat flux should be constant across the interface, which would manifest as a Neumann boundary condition. Currently, I am evaluating the liquid domain, and evaluating the heat flux at the boundary using a forwards second-order finite difference and using that value as the boundary condition for the top of the solid domain. However, this does not preserve the condition that the interface temperature should be constant at $T_{melt}$. My initial thought was to formulate the derivative at the second to last point to match the interface flux, but that method does not seem to be matching the flux as I would expect after computing the temperature. Considering a backwards derivative at the boundary, this is applied as:
$$T_{I-1} = T_{melt} - q^{''}_{interface}\Delta x / -\kappa_l$$
There is also concern regarding the formal order as this approach only uses a first order finite difference approximation and the spatial and temporal order of the CN method should be second order.
I have attempted to solve for the temperature in python using this method and while the temperature looks reasonable, there is still a discontinuity in the heat flux at the interface $\left(\kappa_1 \frac{\partial T_1}{\partial z} \neq \kappa_2 \frac{\partial T_2}{\partial z}\right)$ leading me to believe the application of the boundary equations is incorrect.
