Semi-analytical solution of multilayer diffusion problems with time-varying boundary conditions and general interface conditions

TL;DR

This paper introduces a semi-analytical method for solving complex multilayer diffusion problems with time-varying boundaries and multiple layers, offering improved accuracy and efficiency over existing methods.

Contribution

A novel semi-analytical approach that handles time-varying boundary conditions and multiple layers, surpassing previous methods in accuracy and applicability.

Findings

The method converges faster with more eigenvalues.

It outperforms the unified transform method in accuracy and efficiency.

It can handle problems with many layers and time-varying conditions.

Abstract

We develop a new semi-analytical method for solving multilayer diffusion problems with time-varying external boundary conditions and general internal boundary conditions at the interfaces between adjacent layers. The convergence rate of the semi-analytical method, relative to the number of eigenvalues, is investigated and the effect of varying the interface conditions on the solution behaviour is explored. Numerical experiments demonstrate that solutions can be computed using the new semi-analytical method that are more accurate and more efficient than the unified transform method of Sheils [Appl. Math. Model., 46:450-464, 2017]. Furthermore, unlike classical analytical solutions and the unified transform method, only the new semi-analytical method is able to correctly treat problems with both time-varying external boundary conditions and a large number of layers. The paper is concluded by replicating solutions to several important industrial, environmental and biological applications previously reported in the literature, demonstrating the wide applicability of the work.

Tables

Table 1. Table 1: Relative errors for our semi-analytical solution (Section 3 ) and the unified transform method [ 18 ] (computed using the default value of N x = 15 subscript 𝑁 𝑥 15 N_{x}=15 as per Sheils’ [ 18 ] UTM Heat code) for the test case described in Section 4.3 .

	Runtime	Relative errors
Method	(secs)	$t=0.01$	$t=0.2$	$t=3$
Unified transform	114.93	2.86e-10	4.65e-11	1.98e-11
Semi-analytical [$N=10$]	0.54	7.18e-05	2.26e-06	7.50e-07
Semi-analytical [$N=25$]	1.06	4.56e-06	1.43e-07	1.99e-08
Semi-analytical [$N=50$]	2.43	4.43e-07	1.31e-08	5.86e-09
Semi-analytical [$N=100$]	7.33	5.34e-08	1.70e-09	7.33e-10
Semi-analytical [$N=200$]	15.78	6.55e-09	1.98e-10	9.09e-11
Semi-analytical [$N=300$]	26.66	1.93e-09	5.98e-11	2.70e-11
Semi-analytical [$N=400$]	32.88	8.13e-10	2.56e-11	1.16e-11
Semi-analytical [$N=500$]	44.38	4.15e-10	1.29e-11	6.04e-12
Semi-analytical [$N=600$]	57.02	2.40e-10	7.53e-12	3.21e-12

Equations

$$l_0<l_1<\dots <l_{m-1}<l_m.$$

$$\frac{\partial u_i}{\partial t}=D_i\frac{{\partial }^2{u}_i}{\partial x^2},$$

$$u_i(x,0)=f_i(x),$$

$$a_L{u}_1(l_0,t)-b_L\frac{\partial u_1}{\partial x}(l_0,t)=g_0(t),$$

$$a_R{u}_m(l_m,t)+b_R\frac{\partial u_m}{\partial x}(l_m,t)=g_m(t),$$

$$u_i(l_i,t)=u_{i+1}(l_i,t),$$

$$D_i\frac{\partial u_i}{\partial x}(l_i,t)=D_{i+1}\frac{\partial u_{i+1}}{\partial x}(l_i,t),$$

$$u_i(l_i,t)=u_{i+1}(l_i,t),$$

$${\gamma }_i\frac{\partial u_i}{\partial x}(l_i,t)={\gamma }_{i+1}\frac{\partial u_{i+1}}{\partial x}(l_i,t),$$

$${\gamma }_i\frac{\partial u_i}{\partial x}(l_i,t)=H_i(u_{i+1}(l_i,t)-u_i(l_i,t)),$$

$${\gamma }_i\frac{\partial u_i}{\partial x}(l_i,t)={\gamma }_{i+1}\frac{\partial u_{i+1}}{\partial x}(l_i,t),$$

$$u_i(l_i,t)={\theta }_i{u}_{i+1}(l_i,t),$$

$${\gamma }_i\frac{\partial u_i}{\partial x}(l_i,t)={\gamma }_{i+1}\frac{\partial u_{i+1}}{\partial x}(l_i,t),$$

$${\gamma }_i\frac{\partial u_i}{\partial x}(l_i,t)=H_i({\theta }_i{u}_{i+1}(l_i,t)-u_i(l_i,t)),$$

$${\gamma }_i\frac{\partial u_i}{\partial x}(l_i,t)={\gamma }_{i+1}\frac{\partial u_{i+1}}{\partial x}(l_i,t),$$

$$g_i(t):={\gamma }_i\frac{\partial u_i}{\partial x}(l_i,t)={\gamma }_{i+1}\frac{\partial u_{i+1}}{\partial x}(l_i,t),$$

$$\frac{\partial u_1}{\partial t}=D_1\frac{{\partial }^2{u}_1}{\partial x^2},$$

$$u_1(x,0)=f_1(x),a_L{u}_1(l_0,t)-b_L\frac{\partial u_1}{\partial x}(l_0,t)=g_0(t),{\gamma }_1\frac{\partial u_1}{\partial x}(l_1,t)=g_1(t).$$

$$\frac{\partial u_i}{\partial t}=D_i\frac{{\partial }^2{u}_i}{\partial x^2},$$

$$u_i(x,0)=f_i(x),{\gamma }_i\frac{\partial u_i}{\partial x}(l_{i-1},t)=g_{i-1}(t),{\gamma }_i\frac{\partial u_i}{\partial x}(l_i,t)=g_i(t).$$

$$\frac{\partial u_m}{\partial t}=D_m\frac{{\partial }^2{u}_m}{\partial x^2},$$

$$u_m(x,0)=f_m(x),{\gamma }_m\frac{\partial u_m}{\partial x}(l_{m-1},t)=g_{m-1}(t),a_R{u}_m(l_m,t)+b_R\frac{\partial u_m}{\partial x}(l_m,t)=g_m(t).$$

$$\frac{1}{H_i}{g}_i(t)={\theta }_i{u}_{i+1}(l_i,t)-u_i(l_i,t),$$

$$u_i(x,t)=w_i(x,t)+v_i(x,t),$$

$$a_L{w}_1(l_0,t)-b_L\frac{\partial w_1}{\partial x}(l_0,t)=g_0(t),{\gamma }_1\frac{\partial w_1}{\partial x}(l_1,t)=g_1(t).$$

$$A(t)=\frac{g_0(t)}{a_L}+\left(\frac{b_L}{a_L}-l_0\right)\frac{g_1(t)}{{\gamma }_1},B(t)=\frac{g_1(t)}{{\gamma }_1},$$

$$w_1(x,t)=\underset{{\psi }_{1,1}(x)}{\underbracket{\frac{1}{a_L}}}{g}_0(t)+\underset{{\psi }_{1,2}(x)}{\underbracket{\frac{a_L(x-l_0)+b_L}{{\gamma }_1{a}_L}}}{g}_1(t).$$

$$w_i(x,t)=g_{i-1}(t){\psi }_{i,1}(x)+g_i(t){\psi }_{i,2}(x),$$

$$\frac{\partial v_1}{\partial t}=D_1\frac{{\partial }^2{v}_1}{\partial x^2}+G_1(x,t),$$

$$v_1(x,0)={\tilde{f}}_1(x),a_L{v}_1(l_0,t)-b_L\frac{\partial v_1}{\partial x}(l_0,t)=0,{\gamma }_1\frac{\partial v_1}{\partial x}(l_1,t)=0.$$

$$\frac{\partial v_i}{\partial t}=D_i\frac{{\partial }^2{v}_i}{\partial x^2}+G_i(x,t),$$

$$v_i(x,0)={\tilde{f}}_i(x),{\gamma }_i\frac{\partial v_i}{\partial x}(l_{i-1},t)=0,{\gamma }_i\frac{\partial v_i}{\partial x}(l_i,t)=0.$$

$$\frac{\partial v_m}{\partial t}=D_m\frac{{\partial }^2{v}_m}{\partial x^2}+G_m(x,t),$$

$$v_m(x,0)={\tilde{f}}_m(x),{\gamma }_m\frac{\partial v_m}{\partial x}(l_{m-1},t)=0,a_R{v}_m(l_m,t)+b_R\frac{\partial v_m}{\partial x}(l_m,t)=0.$$

$${\tilde{f}}_i(x)=f_i(x)-\left[g_{i-1}(0){\psi }_{i,1}(x)+g_i(0){\psi }_{i,2}(x)\right],$$

$$G_i(x,t)=D_i\left[g_{i-1}(t){\psi }_{i,1}^{\prime \prime }(x)+g_i(t){\psi }_{i,2}^{\prime \prime }(x)\right]-\left[g_{i-1}^{\prime }(t){\psi }_{i,1}(x)+g_i^{\prime }(t){\psi }_{i,2}(x)\right],$$

$$v_i(x,t)=\sum _{n=0}^{\infty }{c}_{i,n}(t){\hat{\varphi }}_{i,n}(x),$$

Full text

Keywords:   multilayer diffusion; semi-analytical solution; transient boundary conditions; general interface conditions; partition coefficient; jump conditions

1 Introduction

Mathematical models of diffusion in layered materials arise in many industrial, environmental, biological and medical applications, such as heat conduction in composite materials [4, 14, 15], transport of contaminants, chemicals and gases in layered porous media [11, 25], brain tumour growth [1, 13], heat conduction through skin [19], transdermal drug delivery [21, 16] and greenhouse gas emissions [12]. Another important application is in the field of multiscale modelling: for a large number of layers, layered diffusion is one of the simplest examples of a multiscale problem and is ideal for prototyping macroscopic and multiscale modelling approaches [2, 3]. Analytical solutions to such problems are highly valuable as they provide a greater level of insight into the solution behaviour and can be used to benchmark numerical solutions.

The most popular analytical solution approach for multilayer diffusion is classical separation of variables (see, e.g., [14, 7, 23]). By assuming a separated solution in each layer, one immediately finds that the eigenfunctions are coupled via the internal boundary conditions (BCs) at the interfaces between adjacent layers. Together with the external BCs, a system of algebraic equations is obtained, linear in the (unknown) eigenfunction coefficients. With the requirement that this linear system possess non-trivial solutions, one obtains a transcendental equation satisfied by the unknown eigenvalues formulated by setting the determinant of the coefficient matrix of the linear system equal to zero. Since computing the matrix determinant is numerically unstable for large matrices, using this method to compute the analytical solution performs poorly for a large number of layers as either erroneous eigenvalues/roots are returned during the solution procedure or eigenvalues/roots are skipped altogether (as reported by Carr and Turner [2]).

To overcome these issues, Carr and Turner [2] recently developed a semi-analytical solution approach for multilayer diffusion based on the Laplace transform and an appropriately-defined orthogonal eigenfunction expansion. The attractiveness of this approach is that the solution formulas involve a local set of eigenvalues in each layer satisfying simple transcendental equations, resembling those of the single layer problem, which in most cases can be solved explicitly. As a result, the solution performs well for a large number of layers [2] with the approach classified as semi-analytical since computing the inverse Laplace transforms appearing in the solution formulas are carried out numerically.

Recently, Sheils [18] applied the unified transform method, initially proposed by Fokas [5], to the layered diffusion problem and compared the approach to the semi-analytical method of Carr and Turner [2]. Sheils reported that her approach is more accurate, but less efficient and also faulty near the boundaries whenever nonhomogeneous external BCs were applied. We note that her assessment of the accuracy was based on using the default value of 50 terms/eigenvalues in the solution expansions in Carr and Turner [2]’s code. Sheils [18] also remarked that a notable difference between the two methods is that only the unified transform method is applicable to time-dependent external BCs.

In this paper, several key contributions to the literature on multilayer diffusion are presented. Specifically, we:

• (i)

develop a new semi-analytical solution approach for time-dependent external BCs and a general set of internal BCs (Section 3);

• (ii)

study the convergence rate of the new semi-analytical method in (i) and compare it to the classical analytical method (Section 4.2);

• (iii)

carry out a more comprehensive comparison between our semi-analytical method and Sheils’ [18] unified transform method, in terms of accuracy and efficiency, by performing an investigation into the effect of the number of terms/eigenvalues on the solution accuracy (Section 4.3);

• (iii)

explore the effect that changing the interface conditions has on the solution behaviour (Section 4.1);

• (iv)

extend the classical analytical solution approach to a general set of internal BCs by proposing the correct form of the weight function in the orthogonality condition (Appendix C).

The treatment of time-dependent external BCs addresses a deficiency of Carr and Turner [2]’s semi-analytical method, as reported by Sheils [18], and allows for application of the method to a wider range of problems, for example, contaminant transport modelling in layered porous media involving a time-varying inlet concentration [11]. Our method is also not faulty at the end points as is the case for Sheil’s [18] unified transform method whenever nonhomogeneous external BCs are applied. Moreover, general interface conditions permit additional features to be incorporated in the layered diffusion model, such as the volumetric heat capacity in heat conduction problems [8] or the sorption coefficient and partitioning phenomena in chemical transport problems [23].

The remaining sections of this paper are organised in the following manner. Section 2 formulates the multilayer diffusion problem considered in this work while the derivation and implementation of the semi-analytical solution method is outlined in Section 3. Numerical experiments are reported in Section 4. The paper then concludes with a summary and overview of key findings.

2 Multilayer diffusion problem

The one-dimensional multilayer diffusion problem is formulated as follows. Consider a diffusion process on the interval $[l_{0},l_{m}]$, which is partitioned into $m$ subintervals (see Figure 1):

[TABLE]

The resulting domain is represented using the notation $[l_{0},l_{1},\ldots,l_{m-1},l_{m}]$. On each subinterval (layer) $(l_{i-1},l_{i})$, that is for all $i=1,\ldots,m$, we define the diffusion equation:

[TABLE]

where $x\in(l_{i-1},l_{i})$, $u_{i}(x,t)$ is the solution (temperature, concentration, etc.) at position $x$ and time $t$ in the $i$th layer, and $D_{i}>0$ is the diffusion coefficient in the $i$th layer. The initial conditions are given by

[TABLE]

for all $i=1,\ldots,m$, and the external BCs are defined as

[TABLE]

for $t>0$, where the coefficients $a_{L}$, $b_{L}$, $a_{R}$ and $b_{R}$ are non-negative constants satisfying $a_{L}+b_{L}>0$ and $a_{R}+b_{R}>0$, and $g_{0}(t)$ and $g_{m}(t)$ are specified time-dependent functions. Note that the restriction on the coefficients is a classical constraint placed on the single-layer problem ($m=1$) [24], which ensures all of the eigenvalues (see Appendix C) are non-negative and the solution remains bounded as $t\rightarrow\infty$.

Closing the problem requires equations (2)–(5) to be coupled with appropriate internal BCs at the interfaces between adjacent layers: $x=l_{i}$ for all $i=1,\ldots,m-1$. Typically, continuity of the solution and the diffusive flux is implicitly assumed at each interface, that is,

[TABLE]

for $t>0$ and for all $i=1,\ldots,m-1$. In this work, in addition to the above interface conditions (6), we also consider the following choices of interface conditions.

Perfect contact conditions. Consider the following interface condition at the $i$th interface:

[TABLE]

for $t>0$, where $\gamma_{i}>0$. These interface conditions generalise (6) and permit a wider array of problems to be considered as $\gamma_{i}$ and/or $\gamma_{i+1}$ can be different from $D_{i}$ and $D_{i+1}$, respectively. For example, in heat transfer, where $D_{i}$ is the thermal diffusivity $D_{i}:=k_{i}/(\rho_{i}{c_{p}}_{i})$ ($k_{i}$ and $\rho_{i}{c_{p}}_{i}$ are the thermal conductivity and volumetric heat capacity in layer $i$, respectively), the interface conditions (7) allow continuity of the heat flux to be imposed at the interface by setting $\gamma_{i}=k_{i}$.

Jump conditions. Consider the following interface condition at the $i$th interface:

[TABLE]

for $t>0$, where $H_{i}>0$ is the contact transfer coefficient at $x=l_{i}$. These interfaces conditions are a more general form of the perfect contact conditions (7) since dividing equation (8a) by $H_{i}$ yields (7a) in the limit $H_{i}\rightarrow\infty$. For a finite value of $H_{i}$, the interface conditions (8) produce a discontinuity or jump in the solution field at the interfaces. This is useful in applications involving contact resistance at the interfaces, for example in heat transfer, where a very thin resistive (low conductive) layer causes a jump in the temperature [8].

Partition conditions. Consider the following interface condition at the $i$th interface:

[TABLE]

for $t>0$, where $\theta_{i}>0$ is the partition coefficient at $x=l_{i}$. Again, these interface conditions are a more general form of the perfect contact conditions (7), which are given by the special case: $\theta_{i}=1$. The interface conditions (9) maintain a constant ratio between the discontinuous solution values at the interface. This phenomena is common in applications such as analyte transport in porous media [23] and drug release from multi-layer capsules [10].

General interface conditions. Each of the above sets of internal BCs can be expressed in the general form:

[TABLE]

as setting $\gamma_{i}=D_{i}$, $\gamma_{i+1}=D_{i+1}$ and $\theta_{i}=1$, and taking $H_{i}\rightarrow\infty$ produces (6); setting $\theta_{i}=1$ and taking $H_{i}\rightarrow\infty$ yields the perfect contact conditions (7); setting $\theta_{i}=1$ produces the jump conditions (8); and taking $H_{i}\rightarrow\infty$ yields the partition conditions (9).

3 Semi-analytical solution

We now develop our new semi-analytical method for solving the multilayer diffusion problem with time-varying external BCs (2)–(5) subject to the general internal BCs (10).

3.1 Reformulation of problem

Define the interface functions [2]:

[TABLE]

for all $i=1,\ldots,m-1$. Equations (2)–(5) and (10b) can now be reformulated as a sequence of single layer problems [17]:

• •

First Layer ($i=1$):

[TABLE]

• •

Middle Layers ($i=2,\ldots,m-1$):

[TABLE]

• •

Last Layer ($i=m$):

[TABLE]

Clearly, the solution of each of the above problems will involve the interface functions (11), which are unknown. The general idea is therefore to solve the above single layer problems subject to the constraint that the solutions satisfy the, as yet unused, interface condition (10a) [2, 17], which can be rewritten in terms of $g_{i}(t)$ as follows

[TABLE]

for all $i=1,\ldots,m-1$. In summary, solving the single layer problems (12)–(14) subject to the constraint (15) is equivalent to solving the multilayer diffusion problem described by equations (2)–(5) and (10).

3.2 Solution of the single layer problems

To solve each of the single-layer problems (12)–(14), we introduce the substitution:

[TABLE]

for all $i=1,\ldots,m$, where $w_{i}(x,t)$ is chosen so that $v_{i}(x,t)$ satisfies homogeneous versions of the BCs given in equations (12b), (13b) and (14b). For example, in the first layer, $w_{1}(x,t)$ satisfies the BCs:

[TABLE]

If $a_{L}\neq 0$, then one choice222This choice, while not unique, does not effect the uniqueness of the final solution as the initial conditions and source/sink terms are appropriately corrected in (23)–(24). for $w_{1}$ is a linear function in $x$ with time-dependent coefficients. Substituting $w_{1}(x,t)=A(t)+B(t)x$ into (17) and solving for $A(t)$ and $B(t)$ gives

[TABLE]

and hence

[TABLE]

Repeating this process, one sees that, in general, the function $w_{i}(x,t)$ can be expressed as a linear combination of the unknown interface functions at $x=l_{i-1}$ and $x=l_{i}$, that is

[TABLE]

for all $i=1,\ldots,m$. The functions $\psi_{i,1}(x)$ and $\psi_{i,2}(x)$ are identified in equation (18) for $i=1$ and $a_{L}\neq 0$. The remaining cases are summarised in Appendix A. Substituting (16) into (12)–(14), we see that the functions $v_{i}(x,t)$ satisfy the following problems with homogeneous BCs:

• •

First Layer ($i=1$):

[TABLE]

• •

Middle Layers ($i=2,\ldots,m-1$):

[TABLE]

• •

Last Layer ($i=m$):

[TABLE]

The modified initial conditions and the source/sink terms are defined as:

[TABLE]

for all $i=1,\ldots,m$. The solution of each of the single layer problems (20)–(22) can be expressed as an eigenfunction expansion [24]:

[TABLE]

with coefficients and orthonormal eigenfunctions defined as

[TABLE]

The eigenvalues ($\lambda_{i}$) and non-normalized eigenfunctions ($\phi_{i}$) satisfy:

[TABLE]

for $x\in(l_{i-1},l_{i})$ and $i=1,\ldots,m$, subject to the following BCs:

• •

First Layer ($i=1$): $a_{L}\phi_{1}(l_{0})-b_{L}\phi_{1}^{\prime}(l_{0})=0$ and $\phi_{1}^{\prime}(l_{1})=0$.

• •

Middle Layers ($i=2,\ldots,m-1$): $\phi_{i}^{\prime}(l_{i-1})=0$ and $\phi_{i}^{\prime}(l_{i})=0$.

• •

Last Layer ($i=m$): $\phi_{m}^{\prime}(l_{m-1})=0$ and $a_{R}\phi_{m}(l_{m})+b_{R}\phi_{m}^{\prime}(l_{m})=0$.

Note that both the eigenvalues and eigenfunctions are local to each layer and simple to obtain. The form of the eigenvalues $\lambda_{i,n}$ and normalized eigenfunctions $\widehat{\phi}_{i,n}$ for $n=0,1,\ldots$ can be found in [2, Appendix B]. The coefficients $c_{i,n}(t)$ satisfy the ordinary differential equation:

[TABLE]

which is formulated by taking the inner product of both sides of (21a) with $\widehat{\phi}_{i,n}(x)$ and setting:

[TABLE]

The general solution of (25) is given by

[TABLE]

The unknown coefficient $c_{i,n}(0)$ is identified using the initial condition (21b):

[TABLE]

and hence we have that:

[TABLE]

Applying integration by parts to the third and fourth integrals in the expression for $c_{i,n}(t)$ (28) yields:

[TABLE]

where we have set

[TABLE]

We express the above expression for $c_{i,n}(t)$ in terms of the inverse Laplace transform for reasons that will become clear later. Using the convolution property of the Laplace transform allows an equivalent expression to be obtained:

[TABLE]

where, for example, $\overline{g}_{i}(s)=\mathcal{L}\{g_{i}(t)\}$. It follows therefore that our solution approach assumes the Laplace transformations of the boundary functions, namely $\overline{g}_{0}(s)$ and $\overline{g}_{m}(s)$, are known, that is, they exist and can be evaluated analytically or numerically. With $c_{i,n}(t)$ identified, we have solved (20)–(22) for $v_{i}(x,t)$ ($i=1,\ldots,m$). The solution of the single layer problems (12)–(14) are hence given by:

[TABLE]

for all $i=1,\ldots,m$.

3.3 Evaluation of the solution expressions

To evaluate the solution expressions, described by equations (29) and (30), the summations are truncated after a finite number of terms/eigenvalues $N$:

[TABLE]

for all $i=1,\ldots,m$. The inverse Laplace transformations appearing in the coefficients (29) are evaluated using the strategy proposed by Carr and Turner [2], which involves applying the quadrature formula described by Trefethen et al. [22] to the integral representation of the inverse Laplace transform. Given $F(s)=\mathcal{L}\{f(t)\}$, the Laplace transform can be inverted numerically as follows:

[TABLE]

where $\Re\{\cdot\}$ denotes the real part, and $c_{2k-1},z_{2k-1}\in\mathbb{C}$ are the residues and poles of the best $(N_{p},N_{p})$ rational approximation to $e^{z}$ on the negative real line (see, e.g., [22, 2])333Note that because we are dealing with diffusion problems, our interest is in functions that have the functional form $f(t)=\nu_{0}+\sum_{n=1}^{\infty}\nu_{n}e^{-\mu_{n}t}$, where $\nu_{n},\mu_{n}\in\mathbb{R}$. In this case, $F(s)=\nu_{0}/s+\sum_{n=1}^{\infty}\nu_{n}/(s-\mu_{n})$, which yields $F(z_{2k-1}/t)/t=\nu_{0}/z_{2k-1}+\sum_{n=1}^{\infty}\nu_{n}/(z_{2k-1}+\mu_{n}t)$ and therefore the limit of (32) is well behaved as $t\rightarrow 0$ and $t\rightarrow\infty$.. Using (32) gives, for example, the following approximation [2]:

[TABLE]

Note that the quadrature formula (33) requires evaluating $\overline{g}_{i}(s)$, which is the Laplace transformation of the unknown interface function $g_{i}(t)$. To compute these evaluations the solutions are constrained to satisfy the interface condition (15) [2]. Taking Laplace transforms of (15) and rearranging yields:

[TABLE]

where $\overline{u}_{i}(x,s)=\mathcal{L}\{u_{i}(x,t)\}$ is given by:

[TABLE]

The Laplace transformation $\overline{c}_{i,n}(s)=\mathcal{L}\left\{c_{i,n}(t)\right\}$ is linear in the functions $g_{i}(s)$ ($i=1,\ldots,m$)

[TABLE]

With the forms of $\overline{g}_{0}(s)$ and $\overline{g}_{m}(s)$ given, substituting (35) into (34) for all $i=1,\ldots,m$ produces a tridiagonal matrix system in the form:

[TABLE]

where $\mathbf{A}(s)\in\mathbb{C}^{(m-1)\times(m-1)}$, $\mathbf{b}(s)\in\mathbb{C}^{m-1}$ and $\mathbf{x}=\left(\overline{g}_{1}(s),\ldots,\overline{g}_{m-1}(s)\right)^{T}$. The individual entries of $\mathbf{A}(s)$ and $\mathbf{b}(s)$ are given in Appendix B. Solving the linear system (37) evaluated at $s=z_{2k-1}/t$ allows the required evaluations $\overline{g}_{i}(z_{2k-1}/t)$ ($i=1,\ldots,m-1$) to be computed444Note that in the limit as $t\rightarrow 0$ ($s\rightarrow\infty$) and $t\rightarrow\infty$ ($s\rightarrow 0$) the entries of $\mathbf{A}(s)$ and $\mathbf{b}(s)$, which arise from the coefficients of $\overline{g}_{i-1}(s)$ and $\overline{g}_{i}(s)$ and the first term in the expression for $\overline{c}_{i,n}(s)$ (36) (see Appendix B) are well behaved.. For each time $t>0$ at which the semi-analytical solution is sought, the $m-1$ dimensional matrix system (37) must be solved $N_{p}/2$ times.

In contrast to the case of time-independent BCs [2], for $i=0$ and $i=m$, the coefficients (29) feature terms involving inverse Laplace transformations of expressions involving $\overline{g}_{0}(s)=\mathcal{L}\{g_{0}(t)\}$ and $\overline{g}_{m}(s)=\mathcal{L}\{g_{m}(t)\}$. Since the external boundary functions $g_{0}(t)$ and $g_{m}(t)$ are known a priori, these terms can be computed directly using the convolution property, for example:

[TABLE]

For a general and flexible code and to reduce user input, numerically evaluating the integral (38) is preferred. However, preliminary investigation found that the quadrature formula (33) performed better for large $t$ than MATLAB’s in-built integral function, which uses global adaptive Gauss-Kronrod quadrature. Therefore, we use (33) for all $i=0,\ldots,m$.

The only thing left to address is evaluating the unknown interface functions $g_{i}(t)$ ($i=1,\ldots,m-1$) at a given time $t$. Note these evaluations appear in both the expression for $u_{i}(x,t)$ (30) and the coefficients $c_{i,n}(t)$ (29). Since $g_{i}(t)=\mathcal{L}^{-1}\left\{\overline{g}_{i}(s)\right\}$, we can use the approximation (32), giving:

[TABLE]

for $i=1,\ldots,m-1$. Note here that the evaluations of $\overline{g}_{i}(s)$ are the same as those appearing in (33). Finally, we remark that it is due to approximations such as (33) that both the method developed in this paper and the method of Carr and Turner [2] are classified as semi-analytical.

3.4 MATLAB code

The semi-analytical solution developed originally by Carr and Turner [2] is available at

https://github.com/elliotcarr/MultDiff. This repository has now been updated to include the semi-analytical solution developed in this paper, which is applicable to time-dependent external BCs and more general interface conditions.

4 Results and discussion

4.1 Effect of interface conditions

Consider the multilayer diffusion problem described by equations (2)–(5) and (10) with $m=2$ layers, domain $[l_{0},l_{1},l_{2}]=[0,0.5,1]$, diffusivities $D_{1}=1$ and $D_{2}=0.1$, and external boundary condition data $a_{L}=1$, $b_{L}=0$, $g_{0}(t)=1$, $a_{R}=0$, $b_{R}=1$ and $g_{m}(t)=0$. Figures 2a–2d depict the semi-analytical solution over time for four different choices of interface conditions:

• •

Case A: Perfect contact with $\gamma_{1}=D_{1}$ and $\gamma_{2}=D_{2}$.

• •

Case B: Jump conditions with $\gamma_{1}=D_{1}$, $\gamma_{2}=D_{2}$ and $H_{1}=0.5$.

• •

Case C: Partition conditions with $\gamma_{1}=D_{1}$, $\gamma_{2}=D_{2}$ and $\theta_{1}=1.2$.

• •

Case D: Perfect contact with $\gamma_{1}=\gamma_{2}=2.0$.

For Case A, the solution gradient is discontinuous at the interface (see Figure 2a) as $\gamma_{1}\neq\gamma_{2}$ (see interface condition (7b)) while for Case B, the solution is discontinuous at the interface (see Figure 2b) since the transfer coefficient $H_{1}$ is finite. Using equation (8a), the step change in the solution at the interface can be expressed as:

[TABLE]

For given values of $H_{1}$ and $\gamma_{1}$, the difference in solution values at the interface is proportional to the gradient appearing in (40), which explains why the jump discontinuity is absent from the solution at $t=0.01$ (at least visibly), large at $t=0.2$ and small at $t=5.0$ (see Figure 2b). Case C also exhibits a jump discontinuity at the interface with:

[TABLE]

which means that, for a given value of the partition coefficient $\theta_{1}$, the size of the jump discontinuity is directly proportional to the value of $u_{2}(l_{1},t)$. This is confirmed in Figure 2c with the step change in the solution across the interface (41) growing over time. In contrast to the jump conditions of Case B (Figure 2b), the steady state solution is dependent upon the partition coefficient. For Case $D$, both the solution and gradient are continuous at the interface even though the diffusivities $D_{1}$ and $D_{2}$ are not equal (Figure 2d). This is explained by the interface condition (7b) describing continuity of the flux, which reduces to continuity of the gradient if $\gamma_{1}=\gamma_{2}$. Indeed, because of this cancellation, the solutions are invariant under the condition $\gamma_{1}=\gamma_{2}=\gamma$ regardless of the value of $\gamma$.

4.2 Convergence of semi-analytical solution

In this section, we investigate the rate of convergence exhibited by our semi-analytical method by comparing the error against the number of eigenvalues $N$ used in the solution expansions (31). Let $u_{i}(x_{i,j},t)$ denote the exact solution of the multilayer diffusion problem in the $i$th layer evaluated at grid point $x_{i,j}=l_{i-1}+j\Delta x_{i}$ where the grid spacing $\Delta x_{i}:=(l_{i}-l_{i-1})/N_{x}$ and $N_{x}$ is the constant number of divisions in each layer. Furthermore, let $\widehat{u}_{i}^{(N)}(x_{i,j},t)$ denote an approximate analytical solution computed using $N$ terms/eigenvalues in each layer. The relative error of this approximate analytical solution is computed as

[TABLE]

where the maximum is taken over $i=1,\ldots,m$ and $j=1,\ldots,N_{x}+1$. Using the above error definition, we compare the accuracy of the semi-analytical solution (31) to the classical analytical solution derived using separation of variables (see, e.g., [2, 7, 23]). The classical analytical solution subject to the general interface conditions (10) is described briefly in Appendix C.

Figure 3 gives the relative error for both the semi-analytical and analytical solutions versus the number of eigenvalues for Case $C$ of Section 4.1 for both a small and large contrast in diffusivity. For increasing $N$, due to the presence of the exponential $e^{-\lambda_{n}^{2}t}$ in the analytical solution (see Appendix C), the terms in the solution expansion tend to zero extremely rapidly provided $t$ isn’t too small: for $t=0.01,0.2,3$ the exact solution is effectively obtained (i.e., the error falls below the machine epsilon of $2^{-52}$ in MATLAB) after $N=37,8,2$ eigenvalues, respectively.

This behaviour is not observed for the semi-analytical solution and this is the price paid for reformulating the problem (see Section 3.1) to avoid solving a complex transcendental equation arising from a matrix determinant to determine the eigenvalues. While the first term of $c_{i,n}(t)$ (29) approaches zero extremely rapidly for increasing $N$, the second two terms, which arise due to the reformulation, do not. Little difference is observed when comparing Figures 3a and 3b, which demonstrates that the semi-analytical method performs well for both small and large contrasts in the diffusivity.

Both plots in Figure 3 suggest the following linear relationship exists between the logarithm of the error $\varepsilon_{N}(t)$ for the semi-analytical method and the number of eigenvalues $N$:

[TABLE]

or equivalently:

[TABLE]

where $p$ denotes the negative slope of the linear error curves given in Figures 3a and 3b, which is essentially independent across the three values of time $t$. Successive computation of the slopes for increasing $N$ indicates that for each time $t$, the slope $p$ approximately approaches a value of $3$ as $N$ increases. Together these observations suggest the following proposition.

Proposition 1**.**

The convergence rate for the semi-analytical method is

[TABLE]

Proof. This result can be derived by studying the coefficient $c_{i,n}(t)$ appearing in the solution expression (31) and noting that each integral in equation (28) is $O(\lambda_{i,n}^{-2})$ for large $\lambda_{i,n}$. For example, applying integration by parts to the second integral in equation (28) yields:

[TABLE]

for large $\lambda_{i,n}$. Additionally, one can show that each of the constants defined in equation (26) and (27) is either $O(\lambda_{i,n}^{-1})$ or equal to zero. For example, consider the values of $\beta_{i,1,n}$ (26) and $\beta_{i,3,n}$ (27) in the first layer ($i=1$) subject to a Dirichlet BC at $x=l_{0}$. Since:

[TABLE]

(see Appendix A and Appendix B of Carr and Turner [2], respectively), clearly $\beta_{1,3,n}=0$ as $\psi_{1,1}^{\prime\prime}(x)=0$ while:

[TABLE]

for large $\lambda_{i,n}$. Verifying the remaining cases follows similarly. With the above results $c_{i,n}(t)=O(\lambda_{i,n}^{-3})$ for large $\lambda_{i,n}$. Hence, subtracting (31) from (30) gives the following expression for the error:

[TABLE]

for large $\lambda_{i,N}$. By considering the possible cases for the eigenvalues (see Appendix B of Carr and Turner [2]), we observe that $\lambda_{i,n}=O(n)$ for each case involving explicit expressions for the eigenvalues. This is also true for the remaining cases (i.e., the first and last layers under Robin BCs) since the eigenvalues, $\lambda_{1,n}$ and $\lambda_{m,n}$, get closer and closer to $n\pi/(l_{1}-l_{0})$ and $n\pi/(l_{m}-l_{m-1})$, respectively, for large $n$ (large $\lambda_{i,n}$) [20]. It follows that $\lambda_{i,N}^{-3}=O(N^{-3})$ and therefore from (45) we get the desired result. $\blacksquare$

We remark that this slowed convergence is typical of problems with time-dependent BCs. For example, consider the following single-layer problem with a time-dependent Dirichlet BC at $x=0$:

[TABLE]

The solution, truncated after $N$ terms/eigenvalues, is given by

[TABLE]

where the eigenvalues $\lambda_{n}=(2n+1)\pi/2$ and the coefficients are defined as

[TABLE]

For this problem the convergence rate is also $O(N^{-3})$ since $c_{n}(t)=O(\lambda_{n}^{-3})=O(n^{-3})$ due to the presence of the additional integral term involving the boundary function $g_{0}(t)$, which is $O(\lambda_{n}^{-3})$ following (44).

4.3 Comparison to Sheils’ UTM Heat Code

We now compare our MATLAB implementation of the semi-analytical solution derived in Section 3, available at https://github.com/elliotcarr/MultDiff, to the unified transform method of Sheils [18], available at https://github.com/nsheils/UTM_Heat. The chosen test case is an $m=8$ layer version of Case A (from Section 4.1) with domain $[l_{0},l_{m}]=[0,1]$, interfaces $l_{i}=i/m$ for $i=1,\ldots,m-1$, diffusivities $D_{2i-1}=1$ and $D_{2i}=0.1$ for $i=1,\ldots,m/2$, and external boundary condition data $a_{L}=1$, $b_{L}=0$, $g_{0}(t)=1$, $a_{R}=0$, $b_{R}=1$ and $g_{m}(t)=0$. Recall that Sheils’ method is faulty at the end points whenever nonhomogeneous external BCs are applied [18]. To circumvent this issue, we first decompose the solution into its steady state and transient parts: $u_{i}(x,t)=1+v_{i}(x,t)$ for each layer $i$, and solve for $v_{i}(x,t)$, which satisfies homogeneous external BCs.

Table 1 compares the runtimes555All tests cases were carried out in MATLAB R2014b on a MacBook Pro (mid 2014) running MAC OS X Version 10.10.5 with 16 GB of RAM and a 3.0GHz dual-core Intel Core i7 processor. and relative errors of our semi-analytical method (for different numbers of eigenvalues $N$) to Sheils’ unified transform method (with the default solver options). While the semi-analytical method is less accurate for our default number of eigenvalues ($N=50$), as reported by Sheils [18], highly accurate solutions can be computed by simply taking more terms/eigenvalues in the solution expansions (31): with $N=600$ terms, the semi-analytical solution is more accurate for all three times reported and twice as fast as the unified transform method. Moreover, for smaller values of $N$, the semi-analytical solution is very fast with an accuracy that is probably sufficient for most applications.

For time-independent external BCs, one can always homogenise the external BCs before applying the unified transform method (as above), however, it is not clear how to do this for time-dependent external BCs. For example, with $g_{0}(t)=1-e^{-t}$, it is not immediately obvious how to avoid the faulty behaviour of Sheils’ unified transform method at the left end point (Figure 4). This detail, together with the fact that the analytical solution (Appendix C) doesn’t perform well for a large number of layers (as reported by Carr and Turner [2]), leads us to the conclusion that only the semi-analytical method introduced in this paper is able to correctly handle both time-dependent external BCs and a large number of layers.

4.4 Applications

In this section, we apply the semi-analytical solution described in Section 3 to some environmental, industrial and biological applications. The primary aim is to highlight the wide array of problems that can be solved by the code and confirm the validity of the derived semi-analytical solution by reproducing results previously reported in the literature.

4.4.1 Contaminant transport in an aquitard

Liu and Ball [11] define the following governing equations for diffusion of a dissolved contaminant in a layered porous medium:

[TABLE]

for $x\in(l_{i-1},l_{i})$, $t>0$ and $i=1,\ldots,m$, where $C_{i}(x,t)$ is the volume-based aqueous contaminant concentration [$\mu\mathrm{g}/\mathrm{L}$] in the $i$th layer, and $R_{i}$ and $D_{i}$ are constants defined as the dimensionless retardation factor and effective diffusion coefficient in the $i$th layer, respectively. At the interface between adjacent layers ($x=l_{i}$), continuity of mass flux and aqueous concentration is imposed:

[TABLE]

for $t>0$ and $i=1,\ldots,m-1$, where the constant $\epsilon_{i}$ is the porosity in the $i$th layer.

We consider the two-layer test problem described by Liu and Ball [11], where initially $C_{i}(x,0)=0$ in both layers ($i=1,2$). The concentration at the top of the first layer is assumed to be a known but arbitrary function of time, denoted by $f(t)$, while zero mass flux is assumed at the bottom boundary:

[TABLE]

for $t>0$. The function $f(t)$ is assumed to take a Gaussian form:

[TABLE]

where the constant $C_{\mathrm{max}}$ is the peak concentration, and $\mu$ and $\sigma$ are constants. The semi-analytical solution requires the Laplace transformation of $f(t)$, which is found to be:

[TABLE]

Note that the quadrature formula (33) requires evaluation of $\overline{g}_{0}(s)=\overline{f}(s)$ for $s\in\mathbb{C}$. To compute the error function for complex arguments, which isn’t available for the inbuilt MATLAB function erf, we use the erfz function developed by Godfrey [6].

Figure 5 depicts the concentration profiles computed using the semi-analytical solution over time for $C_{\mathrm{max}}=1.0\,\mu\mathrm{g}/\mathrm{L}$, $\mu=2.15\,\mathrm{yr}$, $\sigma=1\,\mathrm{yr}$, $R_{1}=42.42$, $R_{2}=1.67$, $D_{1}=1.6\times 10^{-10}\,\mathrm{m}^{2}/\mathrm{s}$, $D_{2}=2.13\times 10^{-10}\,\mathrm{m}^{2}/\mathrm{s}$ and $\epsilon_{1}=\epsilon_{2}=0.54$. To allow comparison with the figure presented by Liu and Ball [11] the total concentration is shown in Figure 5, which is defined in the $i$th layer as $\widetilde{C}_{i}(x,t)=C_{i}(x,t)\epsilon_{i}R_{i}/{\rho_{b}}_{i}$ [$\mu\mathrm{g}/\mathrm{kg}$], where ${\rho_{b}}_{i}$ is the bulk density of the soil in the $i$th layer with ${\rho_{b}}_{1}={\rho_{b}}_{2}=1.4\,\mathrm{kg}/\mathrm{L}$.

The plot clearly depicts the time-varying inlet concentration, reaching a peak concentration after $2.15\,\mathrm{yrs}$ ($t=\mu$). Note that the total concentration is discontinuous at the interface due to the different retardation factors in both layers. The solutions are in excellent agreement with those given by Liu and Ball [11], giving us confidence that our new semi-analytical method has been formulated correctly for time-dependent external BCs.

4.4.2 Heat conduction in composite materials

Consider the classical heat conduction problem in a composite medium comprising $m$ layers [14], where the temperature distribution in each layer is governed by the heat equation:

[TABLE]

for $x\in(l_{i-1},l_{i})$, $t>0$ and $i=1,\ldots,m$, and interface conditions

[TABLE]

for $t>0$ and $i=1,\ldots,m-1$, where $T_{i}(x,t)$ is the temperature at position $x$ and time $t$ in the $i$th layer and $H_{i}$ is the heat transfer coefficient between layers $i$ and $i+1$. The remaining constants $\rho_{i}$, ${c_{p}}_{i}$ and $k_{i}$ denote the density, specific heat capacity and thermal conductivity, respectively, in the $i$th layer. Note that in general the solution to the above problem cannot be obtained using the semi-analytical method described by Carr and Turner [2] since it is the thermal conductivities $k_{i}$ not the thermal diffusivities $k_{i}/(\rho_{i}{c_{p}}_{i})$ that appear in the second interface condition (54).

As an illustrative example, we consider a classical heat conduction problem that has been solved by several authors, including Mulholland and Cobble [15], Mikhailov et al. [14] and Johnston [9]. The test case considers $m=3$ layers, domain $[l_{0},l_{1},l_{2},l_{3}]=[0,2,4,6]\,\mathrm{cm}$ and the following parameter values: $H_{i}\rightarrow\infty$ for all $i=1,\ldots,3$ (i.e., perfect thermal contact at the interfaces), thermal conductivities $k_{1}=297.64$, $k_{2}=1741.18$ and $k_{3}=565.51$ [units $\mathrm{cal}/(\mathrm{cm}\,^{\circ}\mathrm{C}\,\mathrm{h}$)], densities $\rho_{1}=11.08$, $\rho_{2}=2.71$ and $\rho_{3}=7.4$ [units $\mathrm{g}/\mathrm{cm}^{3}$], and specific heat capacities $c_{p,1}=0.031$, $c_{p,2}=0.181$ and $c_{p,3}=0.054$ [units $\mathrm{cal}/(\mathrm{g}\,^{\circ}\mathrm{C})$]. The initial conditions and external BCs are defined as

[TABLE]

Figure 6, which shows the temperature distribution obtained using the semi-analytical solution at various times, reproduces the results previously reported in the literature [15, 14, 9].

4.4.3 Analyte transport in composite media

Chemical concentration profiles in composite media can be sharply discontinuous at material interfaces due to partitioning phenomena [21, 23]. Trefry and Whyte [23] present the following mathematical model governing analyte transport in a composite medium comprising $m$ laminates/media, consisting of the diffusion equation

[TABLE]

for $x\in(l_{i-1},l_{i})$, $t>0$ and $i=1,\ldots,m$, subject to the following interface conditions:

[TABLE]

for $t>0$ and $i=1,\ldots,m-1$. In the above equations, $C_{i}(x,t)$ is the mobile phase concentration at position $x$ and time $t$ in the $i$th layer, $\alpha_{i}$ is the mobile phase partition coefficient at the interface between layers $i$ and $i+1$ ($i=1,\ldots,m-1$), and $D_{i}$ and $\sigma_{i}$ are the diffusion coefficient for the analyte and the linear sorption coefficient in medium $i$, respectively. Note that the solution of this problem can only be computed using the semi-analytical solution given by Carr and Turner [2] if $\alpha_{i}=1$ for all $i=1,\ldots,m-1$ and $\sigma_{i}=0$ for all $i=1,\ldots,m$.

Consider the test case described by Trefry and Whyte [23] with $m=2$ layers, domain $[l_{0},l_{1},l_{2}]=[0,1,2]$, diffusion coefficients $D_{1}=5$ and $D_{2}=0.05$, partition coefficient $\alpha_{1}=2$ and zero sorption coefficients $\sigma_{1}=\sigma_{2}=0$. Initially, medium $1$ is fully concentrated with analyte, medium 2 has zero initial concentration. A zero flux condition applies at the left boundary and a zero concentration condition applies at the right boundary.

Figure 7 plots the concentration profile for $0\leq t\leq 2$ obtained using the semi-analytical solution. The included three-dimensional plot is in excellent agreement to the one featuring in the paper by Trefry and Whyte [23].

4.4.4 Brain tumour growth

Mantzavinos et al. [13] and Asvestas et al. [1] both consider a one-dimensional reaction-diffusion model for the growth of brain tumours. Due to the heterogeneity of brain tissue, the diffusion coefficient is assumed to be piecewise constant over several regions consisting of either white or grey matter. The spread of malignant cells is governed by the equation [13]:

[TABLE]

for $t>0$ and $i=1,\ldots,m-1$. The equation system (55) is converted into the required form of the multilayer diffusion problem considered in this paper (i.e., without the source term) via the substitution $c_{i}(x,t)=e^{t}u_{i}(x,t)$ [13], which yields the multilayer diffusion problem described by equations (2)–(5) and (6) with $a_{L}=a_{R}=0$, $b_{L}=b_{R}=1$ and $g_{0}(t)=g_{m}(t)=0$.

We consider the test case described by Mantzavinos et al. [13] with $m=3$ regions, domain $[l_{0},l_{1},l_{2},l_{3}]=[-5,-1,1,5]$, diffusion coefficients $D_{1}=D_{3}=0.2$ (grey matter) and $D_{2}=1$ (white matter), and initial point sources of tumour cells at $x=-4$ and $x=2$:

[TABLE]

where $\delta(\cdot)$ is the Dirac delta function. To treat the above initial conditions in our code, we use the approximation:

[TABLE]

with $a=0.1$. Figure 8 plots the cell density over time computed via the semi-analytical solution, replicating the plot presented in [13].

5 Conclusions

This paper has developed a new semi-analytical method for solving multilayer diffusion problems with time-varying external BCs and general internal BCs at the interfaces between adjacent layers. Numerical experiments suggested the semi-analytical solution exhibits a convergence rate of $O(N^{-3})$, where $N$ is the number of eigenvalues used in the solution expansions, a result which was also confirmed analytically. While this is some way away from the exponential convergence rate of the classical analytical solution, the semi-analytical approach possesses other clear advantages such as requiring only simple eigenvalues and performing well for a large number of layers; both due to it avoiding the solution of a complex transcendental equation for the eigenvalues. Numerical experiments demonstrated that solutions can be obtained using the new semi-analytical method that are more accurate and efficient than Sheils’ [18] unified transform method while also not exhibiting faulty behaviour at the external boundaries. Finally, in contrast to classical analytical solutions and the unified transform method, only the semi-analytical method introduced in this paper is able to correctly treat problems with both time-dependent external BCs and a large number of layers.

Acknowledgments

EJC acknowledges funding from the Australian Research Council (DE150101137).

Appendix A Psi functions

The functions $\psi_{i,1}(x)$ and $\psi_{i,2}(x)$ described in equation (19) are defined as follows:

• •

First Layer ($i=1$)

[TABLE]

• •

Middle Layers ($i=2,\ldots,m-1$)

[TABLE]

• •

Last Layer ($i=m$)

[TABLE]

Appendix B Linear system

This appendix formulates the entries of the matrix $\mathbf{A}(s)$ and $\mathbf{b}(s)$ featuring in the tridiagonal matrix system (37). Let $a_{i,j}$ denote the $(i,j)$ entry of $\mathbf{A}$ and $b_{i}$ denote the $i$th entry of $\mathbf{b}$.

Consider equation (36) and define $c_{i,n}^{(k)}(s)$ ($k=1,2,3$) as follows:

[TABLE]

Substituting (35) into (34), and rearranging identifies

• •

the subdiagonal of $\mathbf{A}$:

[TABLE]

• •

the diagonal of $\mathbf{A}$:

[TABLE]

• •

the superdiagonal of $\mathbf{A}$:

[TABLE]

The entries of $\mathbf{b}$ are given by:

[TABLE]

Appendix C Analytical solution

Using separation of variables, the following analytical solution of the multilayer diffusion problem described by equations (2)–(5) and (10) can be derived [23, 2, 7] when $g_{0}$ and $g_{m}$ in equations (4)–(5) are independent of time $t$:

[TABLE]

where $w_{i}(x)$ denotes the steady-state solution in the $i$th layer (see, e.g., Carr and Turner [2]). The eigenvalues ($\lambda$), which form a global set valid across all layers, and the non-normalized eigenfunctions ($\phi_{i}$) satisfy a series of coupled Sturm Liouville problems involving homogeneous versions of the internal and external BCs:

[TABLE]

Substituting the form of the eigenfunctions:

[TABLE]

into the above internal and external BCs yields a linear system, which can be expressed in matrix form as $\mathbf{A}(\lambda)\mathbf{x}=\mathbf{0}$, where $\mathbf{x}=(\zeta_{1},\xi_{1},\ldots,\zeta_{m},\xi_{m})^{T}$ and $\mathbf{A}\in\mathbb{R}^{2m\times 2m}$. The notation $\zeta_{i}(\lambda)$, $\xi_{i}(\lambda)$ and $\mathbf{A}(\lambda)$ is used since they each depend on $\lambda$. The eigenvalues $\lambda_{n}$ ($n=0,1,\ldots$) are the non-negative roots of the transcendental equation:

[TABLE]

For each eigenvalue $\lambda_{n}$ ($n=0,1,\ldots$), an eigenfunction $\phi_{i,n}(x):=\phi_{i}(\lambda_{n};x)$ is defined, where the coefficients $\zeta_{i}(\lambda_{n})$ and $\xi_{i}(\lambda_{n})$ ($i=1,\ldots,m$) are determined by finding a non-trivial solution of $\mathbf{A}(\lambda_{n})\mathbf{x}=\mathbf{0}$ [2]. The eigenfunctions are orthogonal over the full domain $[l_{0},l_{m}]$ with respect to the weight function $p_{i}(x)=\gamma_{i}D_{i}^{-1}\prod_{k=1}^{i-1}\theta_{k}$. A proof of this result follows closely the proof given by Trefry and Whyte [23] for $\gamma_{i}=D_{i}$. Hence, using the initial condition, we have that:

[TABLE]

where $\tilde{f}_{i}(x)=f_{i}(x)-w_{i}(x)$.