Transient Diffusion in Bi-Layer Composites With Mass Transfer Resistance: Exact Solution and Time Lag Analysis

Exact analytical and closed-form solutions to the transient diffusion in bi-layer composites with external mass transfer resistance are reported. Expressions for the concentrations and the mass permeated are derived in both the Laplace and time domains through the use of the Laplace transform Inversion Theorem. The lead and lag times, which are often of importance in the characterization of membranes and arise from the analysis of the asymptotic behavior of the mass permeated through the bi-layer composite, were also derived. The presented solutions are also compared to previously derived limiting cases of the diffusion in a bi-layer with an impermeable wall and constant concentrations at the upstream and downstream boundaries. Analysis of the time lag shows that this membrane property is independent of the direction of flow. Finally, an outline is provided of how these transient solutions in response to a step function increase in concentration can be used to derive more complex input conditions. The importance of adequately handling boundary layer effects has a wide array of applications such as the study of bi-layers undergoing phenomena of heat convection, gas film resistance, and absorption/desorption.


INTRODUCTION
The mathematical description of the transient behavior of a composite membrane and, in particular, bi-layer composite membranes, has important implications in drug release devices from planar matrix devices (Cabrera et al., 2006;Cabrera and Grau, 2007), in drug delivery through the skin (Ghanem et al., 1987;Couto et al., 2014), in the study of the cornea (Cooper and Kasting, 1987;Couto et al., 2014), studying the role of membranes in keeping compounds from contaminating the surroundings (Kalbe et al., 2002;Edil, 2003), in vacuum insulation panels that creates high performance thermal insulation (Garnier et al., 2011), mass diffusion of neutral species in the fabrication of multi-layer thin-film (Goldner et al., 1992), metallic thermal protection system (Gu et al., 2016), semiconductor composite in bi-layer organic solar cells (Hatton et al., 2007), biodegradable bi-layer film barrier using gelatin and chitosan (Rivero et al., 2009), gas transport in inorganic/organic hybrid structures (Jang and Han, 2009), and gas barrier of organic-inorganic hybrid coatings (Minelli et al., 2010), When modeling the phenomena of heat conduction or gas permeation, the exchange of heat or chemical at the surfaces contacting with the medium are often assumed to occur "infinitely" fast. Mathematically, this translates into the assumption that the temperature (for heat conduction) or concentration (for mass transfer) of the bi-layer at the surface equates that of the contacting medium, simplifying greatly the calculations necessary to find an analytical solution. The contacting medium is generally considered to be maintained at a constant temperature or concentration. To account for boundary effects as shown in Figure 1A, whether this means heat convection, diffusion-limited reaction at the surface or mass transfer resistance, the constant boundary condition can be altered to: with C ∞ being the concentration of the solution and assuming that the volume of the contacting solution is of large enough volume that its concentration is invariant. For example, in fuel resistance of rubber where phenomena of absorption of the fluid on the exposed side of the rubber membrane and desorption on the other side play important roles the transient behavior of this membrane (Forte and Leblanc, 1992). Frisch showed the benefits of a boundary-layer resistance in the form of a protective film layer on top of polymer membranes (Frisch et al., 1984). Boundary effects are seen in food packaging with polymer films to reduce the permeation of oxygen or water (Smith and Peppa, 1991;Laoubi and Vergnaud, 1996). A film barrier can also form during dispersed-drug release from planar matrices (Zhou and Wu, 2002) that, without adequate accounting of these boundary effects, leads to an underestimation of the drug amount released into the medium. Another important implication of boundary layers appear in the context of rotating disc electrodes (Levich, 1962;Gough and Leypoldt, 1980).
There have been considerable efforts in developing solutions to bi-layer systems involving the mass-transfer limited boundary conditions (Equation 1). Ramkrishna and Amundson (Ramkrishna and Amundson, 1974), as well as (Mikhailov et al., 1983), have proposed solution schemes based on solving the eigenvalues and eigenfunctions of a Sturm-Louiville type system to study the transient heat conduction behavior of composites undergoing convection type boundary conditions. Antonopoulos and Tzivanidis (Antonopoulos and Tzivanidis, 1996) outlined a procedure to evaluate transient heat conduction in n-layer composites with a convection type of boundary condition combining the separation of variables technique with an orthogonal expansion of the functions. Yet, the analytical, as opposed to numerical, evaluation of the eigenvalues, eigenfunctions, and the integral relationships arising from these approaches generally remains challenging for a bi-layer composite undergoing convective effects. This technique was further developed to investigate the approximate behavior of heat conduction with combined radiative and convective boundary conditions (Miller and Weaver, 2003) through a linearization of the radiative effects. An 'analytical' method of this problem was also developed by de Monte (de Monte, 2000;de Monte, 2002) based on Vodicka's approach and applied to one-dimensional media to study unsteady heat conduction and extending to n-layers. We also note related works by Goldner et al. (1992) who obtained the analytical equations to a bi-layer diffusion system subjected to an impenetrable wall boundary using the Inversion theorem of Laplace transforms. Subramanian and White (Subramanian and White, 2001) derived an exact solution to a bi-layer composite undergoing galvanostatic boundary conditions using a modified separation of variables method.
In addition to the importance of modeling the transient response in a permeation process, especially for slowly diffusing processes, the "integral permeation" method has played a key role in the determination of the physical properties of a membrane, and in particular, its constitutive diffusion coefficients and permeability constants (Rutherford and Do, 1997). By dividing the permeation process into its transient and steady-state components as shown in Figure 1B, one can represent the transient component by a "time lag" parameter. This concept was introduced in the work of Daynes (1920) and Barrer and Rideal (1939). These findings outlined the now widely used relationship of the time lag (t lag ) for a single layer with assumptions of constant concentrations maintained on either side of a single-layer material: with L being the thickness and D being the diffusion coefficient of the material. Also in Figure 1B, the steady-state components of the upstream and downstream mass permeated are extended using the dashed lines and their x-axis intercepts represent the lead time t + and the time lag t L . The forward mean first passage time is provided by t + t L − t + . The time lag technique has since been expanded to permeation with non Fickian diffusion, polymer membranes, reaction, metallic solids, through catalyst particles among other things (Rutherford and Do, 1997). Notably, there has been an extended amount of work on the characterization of composite membranes using the time lag approach (Jaeger, 1950;Frisch, 1957Frisch, , 1958Frisch, 1959;Pollack and Frisch, 1959;Ash et al., 1963Ash et al., , 1965Siegel, 1991;Liang and Siegel, 2006;Al-Qasas et al., 2014). In this work, the closed-form and exact transient solutions to a single-layer material and a bi-layer composite undergoing mass transfer resistance at its outside boundaries are derived using the Laplace Inversion Theorem (Jaeger, 1950). The time lag and steady-state components of the transient response to an input step function are also derived. This work aims to support the characterization of bi-layer composites that have significant boundary layer effects. The presented solutions are also shown to be a generalization to various existing solutions in the literature. In addition to the similarities between diffusion and heat conduction, we also show how the solution forms presented in this work can be used to study the transient responses of bi-layer composites to inputs more complex than the usual step function.

THEORETICAL ANALYSIS
A schematic of the problem is shown in Figure 1A. For both the single-layer and bi-layer derivation, the flow within the composite is assumed to follow Fick's second law with mass transfer resistance occurring at the external boundaries. The assumed direction of the flow is from left to right.

Diffusion in a Single-Layer with Transfer Resistance
Let us note that the following derivation of the single-layer is equivalent to that of Ash (2001). The key differences with that work are the use of different constant definitions and the application of the nondimensionalization of the parameters. This subsection was written in an effort to provide results that are consistent between the previously derived single-layer and the newly derived bi-layer solutions as well as making comparisons with other results in the literature under one set of notations.

Single-Layer Diffusion with Upstream and Downstream Transfer Resistance
The phenomenon of diffusion is modeled using "Fick's 2 nd Law," i.e., Crank (1979): subject to the boundary and initial conditions: and The convention C(x,t) is used throughout the paper. B 1 and B 2 are the Biot numbers for systems with heat conduction undergoing heat convection at the outside boundaries, the mass transfer Biot numbers in mass diffusion processes or the reaction constant for a first-order reaction. The Biot number is a quantity that relates the inner diffusion to the mass transfer resistance at the surface. Let's note that when both B 1 and B 2 → ∞, these conditions reduce the problem to that with the analytical solution derived by Daynes (1920). Dimensionless variables for the concentration, time, and position are introduced: The dimensionless quantity τ is also known as the mass Fourier number, which can be thought as the ratio between the diffusive transport rate and the storage rate. By using the introduced dimensionless variables and transforming Equation 3 into the Laplace domain, the relation with time is transformed into a dependence on the complex variable s: The solution for c, the Laplace-domain concentration, assumes the form (q s √ ): c a 1 e qy + a 2 e −qy (10) and is subjected to the transformed boundary and initial conditions: and d c dy Solving for c yields the exact Laplace domain solution: The mass permeated can be calculated as follow: and the corresponding dimensionless mass permeated is provided by: The dimensionless mass permeated expressions for the upstream and downstream fluxes, and Applying the Inversion Theorem for Laplace transforms (Jaeger and Carslaw, 1959), we find that the eigenvalues of the timedomain solution are given by the roots of: λ n (B 2 cos(λ n ) − λ n sin(λ n )) + B 1 (B 2 sin(λ n ) + λ n cos(λ n )) 0 and the time-domain solutions provided by Equations 19-21.
The time lag as described by Siegel (1991) for the single-layer problem was given by: Using the work in Siegel (1991), the lead time was also derived:

Single-Layer Diffusion With Constant Upstream and Downstream Concentrations
By letting B 1 → ∞ and B 2 → ∞, the boundary conditions simplify to c(0) 1 and c(1) 0. We show that one can rederive the expressions found by Daynes (1920) with the equations for the eigenvalues: concentration: and downstream mass permeated: Noting that letting τ → ∞ and m down 0, one gets that τ 1 6 or, equivalently, that t L 2 6D .
These expressions are equivalent to the solutions presented in Gough and Leypoldt (1980). The governing equations of the bi-layer system shown in Figure 1A are provided by the Fick's 2 nd law: subject to the external boundary conditions: and the interface boundary conditions: The upstream interface is set at x 1 L 1 and the downstream interface is located at x 2 L 2 . Equation 32 does not contain a second term on the right-hand side (assumes that C ∞ 0) due to C 0 accounting for the difference between the upstream (left) and downstream (right) side concentrations of the contacting medium, as only the gradient makes up the driving force of the permeation process. One can then simply shift the solutions by the value of the downstream medium concentration. The interface boundary condition in Equation 34 also allows for a mismatch between the concentrations between the first and second layers of the bilayer through the constant parameter K. The following nondimensional variables and parameters are introduced: Noting that q s √ , with the complex variable s in the Laplace domain, equivalent expressions are derived: with the transformed boundary conditions: The constant ρ 2 results from the nondimensionalization of Fick's second law for the second layer using the dimensionless time introduced in the nondimensionalization of the first layer. This constant is also the ratio between the mass Fourier number of the first layer and that of the second layer. κ arises from the nondimensionalization of the flux at the interface between the first and second layer. It relates to the ratio between the Biot numbers of the two layers with the assumption that the convective term is the same for both layers. Thus, ρ and κ are dimensionless quantities that quantify how the two layers relate to one another. The thickness of each layer is normalized such that the dimensionless thickness is 1. By solving the equations for the variables c 1 and c 2 , the exact solutions in the Laplace domain are obtained by Equations 43 and 44.
with the definitions: and respectively. Applying the Inversion theorem, the eigenvalue relationship is given by Equation 49  sin(λ n ) cos ρλ n B 1 B 2 κρ − ρλ 2 n + cos(λ n ) sin ρλ n B 1 B 2 − κρ 2 λ 2 n − sin(λ n ) sin ρλ n B 1 κρ 2 λ n + B 2 λ n + cos(λ n ) cos ρλ n B 1 ρλ n + B 2 κρλ n 0 Frontiers in Chemical Engineering | www.frontiersin.org January 2021 | Volume 2 | Article 605197 2B 1 λ n cos y 1 λ n B 2 sin ρλ n + ρλ n cos ρλ n −κρ sin y 1 λ n B 2 cos ρλ n − ρλ n sin ρλ n e −λ 2 n τ D (50) 2B 1 λ 2 n −sin(λ n ) B 2 sin ρλ n +ρλ n cos ρλ n +κρcos(λ n ) B 2 cos ρλ n +ρλ n sin ρλ n 1−e −λ 2 D sin(λ n ) sin ρλ n B 1 B 2 κρ 2 + 1 + B 1 κρ 2 + B 2 − ρ 2 λ 2 n (1 + κ) −ρ cos(λ n ) cos ρλ n B 1 B 2 (κ + 1) + B 1 + B 2 κ − λ 2 n 1 + κρ 2 +ρλ n sin(λ n ) cos ρλ n B 1 κρ 2 + 1 + B 2 (1 + κ) + 2 +λ n sin ρλ n cos(λ n ) B 1 ρ 2 (κ + 1) + B 2 1 + κρ 2 + 2κρ 2 (54) Two representative examples of the transient concentration expressions, c 1 and c 2 , are presented in Figures 2A,B. The first notable distinction between the two examples is that (a) settles more quickly to its steady-state value. The primary reason for this is that (b) has a stronger mass transfer resistance occurring at the downstream boundary, acting almost like an impermeable wall that limits the flux through the bi-layer. Additionally, in (b), due to κ being of value 1, there is no sharp difference in concentration behavior at the interface between the two layers. In the case of (a), the value of 0.1 for κ leads to a sharp contrast in the slope of the concentrations between the two layers. Finally, the major contrast between accounting for mass transfer resistance or not, lays in the fact that the dimensionless concentrations at the upstream boundary is not 1 nor is it 0 at the downstream boundary when the values of B1 and B2 are finite values.

Bi-Layer Diffusion With An Impermeable Wall at One End
There are a number of exact solutions that were previously derived for both the single-layer and bi-layer materials.
This solution not only assumes no flux at the downstream boundary but also that the transfer at the upstream boundary happens so quickly that the concentration there is essentially constant. We note here that, when compared with the previous work on this solution by Goldner et al. (1992), the eigenvalue expression (Equation 57) looks slightly different. In the attempt by the authors to simplify Equation 57 and divide the eigenvalue expression by cos(λ n )cos(ρλ n ), the authors omitted the possibility that, in some cases, the eigenvalues can lead to term containing the cosines to be 0, and thus leading to the tangent to be +∞. However, this is only problematic in the scenario that such eigenvalue would also make sin(λ n ) or sin(ρλ n ) to be 0. A representative example of this type of bi-layer is shown in Figure 2C. Most notably, the concentration is uniform across the bi-layer at steady-state. If B 1 is a finite value, then this steady-state concentration would be less than 1. The κ value of 5 here makes it such that the slope of the concentration in the first layer is sharper than the second layer. This difference, however, fades over time as the solution converges to steady-state.
An example of this well-known case is shown in Figure 2D. The main characteristics are that the upstream concentration is always 1, the downstream concentration is always 0.

Extension of the Step Response Solutions to Other Types of Input Responses
Using the Laplace expressions for both the single-layer and bilayer systems, the derivation of other types of input (e.g., oscillatory, exponential decay, etc.) can be found trivially by factoring out 1 s from the expressions in Equations 19,20,43, and 44 and replacing it with the Laplace s-domain function representing the input of choice. Numerically inverting Laplace solutions into the time-domain can be done through various methods, some of which are outlined in Abate and Whitt (2006). Alternatively, one can apply the truncated time domain convolution: While the time-domain expressions may be convoluted in nature, the time-dependent elements of the solutions are relatively simple. Thus, the convolution can be straightforwardly calculated for many cases such as sine, cosine, or exponential inputs.
Two representative examples are shown in Figures 2E,F. In the first case (e) with u(t) e −0.1τ , the exponential decay function starts at a high value, leading to an increase in concentration throughout the bi-layer. As the driving force of the input decays over time toward 0, the concentrations fall from high values around τ 1.69 and eventually approaches 0 across the two layers. In the second case scenario with u(t) 1 + cos(2τ), the behavior of the bi-alyer for a more complex input made up of a step function component and an oscillatory component is shown. This shows that the step function solutions can be easily adapted to handle more complex behaviors of the upstream "input" to the system.

Time Lag Analysis
The steady-state portion of m can be obtained from the steadystate portion of its transform, lim s → 0 m(s) as shown in Jaeger (1950). However, in this case, it remains quite challenging to derive for the bi-layer problem using that technique. Instead, the method outlined in Siegel (1991) is applied that makes analogies to the concepts in electrical engineering treating layers in a composite material like resistances in series. By assuming that the upstream and downstream boundaries are treated as layers of no thickness, the permeability P of the first and second layers take the following forms: with (01) representing the first layer and (23) representing the second layer. Using previous results that were obtained in the work of Siegel (1991), the time lags of the uncombined two layers are provided by: Let's note that the time lags are nondimensionalized such that τ D1t L 2 1 . The following expressions of forward mean first passage time of the second layer τ (23) + (defined as τ + τ L − τ + , with τ + being the lead time) as well as the backward mean first passage time of the first layer τ (01) − (defined as τ − τ L − τ − , with τ − being the backward lead time) are also introduced to derive the time lag of the bi-layer. Applying the time lag (τ L ) expression for a series composite Siegel (1991), the following expression for τ L is obtained: Similarly, the relationship in Siegel (1991) for the mean passage time in a bi-layer composite is provided by: and τ + 1 2 + ρ 2 1 2 The lead time τ + is then simply given by: (73) Figure 3 shows the time lags for both the single-layer and bi-layer cases. In the single-layer case (a), which only depends on dimensionless parameters B 1 and B 2 , an axis of symmetry (shown as a dotted white line) can be drawn across the heatmap. It indicates that the time lag for the single-layer system is independent as to whether the flow from the left to right or right to left. A similar observation to the invariance of the time lag with respect to flow direction can be made for the bi-layer system and is shown in (c) and (d). However, since the quantities κ, ρ, and τ L are all dependent on how the layers are ordered. Given that ρ and κ are simply ratios between the first and second layer, one simply needs to take the reciprocal of these values to reverse the flow direction. For τ L , one needs to transform from the Fourier number of the first layer to that of the second layer, which can be simply done through dividing τ L by ρ 2 or τ ′ L as in Figure 3D.
In addition to the invariance to flow direction, it is also interesting to point out that the time lag is greatly affected by the values of B 1 and B 2 . If both of these values are low, the time lag can be extremely large which coincides with the expected behavior that having large mass transfer resistance at the surfaces is not good for reaching steady-state permeation. This also shows that, depending on the values of κ and ρ in the bi-layer case, it may be easier to bring down the time lag by increasing one of the Biot number rather than the other. This can be seen in Figure 3C, for which it is easier to bring the time-lag down by increasing B 2 . At the other end of the spectrum, high values for both B 1 and B 2 yield low time lag values.
In Figure 4, the steady-state and transient solutions for the dimensionless mass permeation expressions are combined for  . Evidently, the slope of the upstream and downstream steady-state solutions are the same since the flux coming into the membrane needs to exit at an equal rate or there would be accumulation within the bi-layer. This type of graph is important because, oftentimes, it is easier in practice to measure the mass permeated rather than the concentration value at a fixed location. Given time series data of the mass permeated, one can use the steady-state and transient solutions provided in this work to calculate missing membrane property values of the materials. In presented examples, it can be seen that the solution that reaches steady-state has high values of both B 1 and B 2 , an observation that can also be seen in Figure 3. One can also use the slopes of the steady-state solutions to find missing values by noting, for example, that the steady-state slopes of the mass permeated expressions do not depend on κ and proceed through process of elimination.

CONCLUSION
There has been a large body of literature treating the problem of diffusion in n-layer composites. Yet, the derivation of closed-form solutions remain seldom and often limited to single-layer or bilayer systems due to the difficult and tedious nature of finding such solutions. In this work, we presented analytical solutions to the phenomena of diffusion in a bi-layer with mass transfer resistance at the external boundaries for each of the layers as well as the corresponding expressions for the asymptotes of the upstream and downstream mass permeated. We showed how these relates to previously derived solutions of diffusion in a bilayer with an impermeable wall as well as diffusion in a bi-layer with constant concentrations on either side. The step response solutions are also extended to time-varying inputs such as exponentially decaying and oscillating inputs. Notably, both the time-domain and Laplace domain solutions to the bi-layer can be used through numerical inversion of the Laplace solution or by applying a convolution on the provided time-domain solutions.
Also, the mean first passage time that is the difference between the lead and lag times does not depend on B 1 . Since the lead time is dependent on the time lag, however, it is also dependent on all the four dimensionless quantities B 1 , B 2 , κ, and ρ. It was also shown that flipping the membrane materials around of the system has no bearing on the time lag for the single-layer and the bi-layer.
When comparing the generalized solutions with the simplified equations with constant concentrations on either side of the bilayer membrane, there are notably two additional parameters in B 1 and B 2 . It is likely that, in the absence of prior knowledge about the system, one cannot easily fit all the parameters at once from experimental data, especially when accounting for measurement errors. Knowledge about B 1 and B 2 can be derived most easily from studying the transient mass permeation in single-layer experiments of each of the components of the bi-layer. This is possible because the nondimensionalized lag time and the slope of the mass permeation in a single-layer membrane only depend on B 1 and B 2 . Using experimental data from the bi-layer, one can, for example, use the fact that the slope of the mass permeation expressions are both independent of ρ. In this scenario, one can then use the lag time τ L of the bi-layer to find κ. Real quantities of thickness, L 1 and L 2 , are usually available or can be measured. Finally, the transient behavior that is provided by the derived equations of the generalized solutions to the bi-layer can be used to fit missing parameters. The objective of this work is to bring more insights as to how to characterize the permeation through bi-layers having nonnegligible mass transfer resistances at its external boundaries by providing both the transient and steady-state analytical solutions.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
JM conceptualized the initial idea of the paper. AT expanded these notes into a full manuscript and derived the equations presented here. ES provided help in the writing and support toward the completion of this work. All authors contributed to manuscript revision, read, and approved the submitted version.