How do I set up the tridiagonal matrix for a heat diffusion with layers of different thermal diffusivity?

Question

I have Scala code that recreates the Crank-Nicolson solutions for the diffusion equations, and matches 'Excel for Scientists and Engineers' (Joe Billo, Wiley).

However, I would like to be able to replicate the results from H. Asan's 'Numerical computation of time lags and decrement factors' (including the decrement factor with very thin materials) and also the multi-layer positioning results from the BS2013 paper 'Optimizing insulation-thermal mass layer distribution from maximum time lag and minimum decrement factor point of view'.

Both papers suggest that Crank-Nicolson is at the heart of the approach, however there is a complexity with the treatment of energy flow at the outer and inner surface, and with different material layers, that are not covered in the simple diffusion treatments that I found.

How do I set up the matrix - or how do I connect a series of such matrices together?

I had thought that perhaps a simple treatment would give a given wall thickness a behavior like a delay factor and an amplitude scaling - but it does not seem to be possible to do this and sum the phase delays and multiply the decrements, because the results above suggest that the ordering is not transitive.

How can I proceed - given that I'm a lot better at coding than solving PDEs.

Answer 1

The heat equation for non-constant coefficients, which is more or less what you're doing here, takes the form, $$ \frac{\partial\psi}{\partial t}=\nabla\cdot\left(\kappa\nabla\psi\right) $$ We can simplify by using a single dimension: $$ \frac{\partial\psi}{\partial t}=\frac{\partial}{\partial x}\left(\kappa\frac{\partial\psi}{\partial x}\right) $$ Rather than applying the derivative operator to the term in the parenthesis (to get a diffusive term & an advective term), we discretize here (though I heartily encourage you to work through the math to get this next equation) to obtain \begin{align} \frac{\psi^{n+1}-\psi^n}{dt}&=\frac1{2dx^2}\left(\kappa_{i+1/2}\left[\psi_{i+1}^{n+1}-\psi_i^{n+1}\right]-\kappa_{i-1/2}\left[\psi^{n+1}_i-\psi^{n+1}_{i-1}\right]\right)\\ &+\frac1{2dx^2}\left(\kappa_{i+1/2}\left[\psi_{i+1}^{n}-\psi_i^{n}\right]-\kappa_{i-1/2}\left[\psi^{n}_i-\psi^{n}_{i-1}\right]\right) \end{align} where $$ \kappa_{i+1/2}=\frac{1}{2}\left(\kappa_{i+1}+\kappa_i\right) $$ You can solve for the $\psi^{n+1}$ terms, $$ -\delta\kappa_{i+1/2}\psi^{n+1}_{i+1}+ \left(1+\delta\left(\kappa_{i+1/2}+\kappa_{i-1/2}\right)\right)\psi^{n+1}_i-\delta\kappa_{i-1/2}\psi^{n+1}_{i-1}=RHS $$ where $\delta=dt/2dx^2$ and $RHS$ contains the $\psi^n$ terms; this is solved as the normal constant-coefficient case. However, you would need to generate the $\kappa(x)\to\kappa_i$ mapping, just as you do for $\psi(x)\to\psi_i$.

You can test this by letting \begin{align} \kappa(x)&=1+\left(\frac x5\right)^2 \\ \psi(x,0)&=\begin{cases}15 & x<0 \\ 20 & x=0 \\ 25 & x>0\end{cases} \end{align} As you let $\psi(x,t)$ evolve, you'll find that the solution tends to the stationary form $$ \psi(x,t)=20\left(1+\frac1\pi\tan^{-1}\left(\frac x5\right)\right) $$ I ran this to $t=4$ and got decent results, you may want to run longer. The above solution can be derived by letting $\partial_t\psi=0$ and solving the resulting ODE by normal techniques.