Generalized Laplace-Type Transform Method for Solving Multilayer Diffusion Problems

Multilayer diffusion problems have found significant importance that they arise in many medical, environmental, and industrial applications of heat and mass transfer. In this article, we study the solvability of a one-dimensional nonhomogeneous multilayer diffusion problem. A new generalized Laplace-type integral transform is used, namely, the Mρ,m -transform. First, we reduce the nonhomogeneous multilayer diffusion problem into a sequence of one-layer diffusion problems including time-varying given functions, followed by solving a general nonhomogeneous one-layer diffusion problem via the Mρ,m-transform. Hence, by means of general interface conditions, a renewal equations’ system is determined. Finally, the Mρ,m-transform and its analytic inverse are used to obtain an explicit solution to the renewal equations’ system. Our results are of general attractiveness and comprise a number of previous works as special cases.

As epidemiological models, reaction-diffusion problems are widely used to model and analyze the spread of diseases such as the global COVID-19 pandemic caused by SARS-CoV-2. These models describe the spatiotemporal prevalence of the viral pandemic and apprehend the dynamics depending on human habits and geographical features. The models estimate a qualitative harmony between the simulated prediction of the local spatiotemporal spread of a pan-demic and the epidemiological collected datum (see [22,23]). These data-driven emulations can essentially inform the respective authorities to purpose efficient pandemicarresting measures and foresee the geographical distribution of vital medical resources. Moreover, such studies explore alternate scenarios for the repose of lockdown restrictions based on the local inhabitance densities and the qualitative dynamics of the infection. For more applications, one can refer, e.g., to [24,25].
Although the numerical methods are usually applied to solve the diffusion problems, especially in the heterogeneous permeable media, the analytic solutions, when available, are characterized by their exactness and continuity in space and time. In the context of obtaining numerical solutions for such models, we refer to the following references [26][27][28][29][30]. In this work, we focus on analytic solutions of certain nonhomogeneous diffusion problems in multilayer permeable media. Here, the retardation factors are assumed to be constant, the dispersion coefficients vary across layers, but being constants within each layer, and the free terms are (arbitrary) time-varying functions.
Analytic and semianalytic solutions of multilayer diffusion problems are developed by using the integral transforms [6,25,[31][32][33][34][35][36][37][38][39][40]. Applying Laplace transforms, to solve multilayer diffusion problems, has advantages as an applicable tool in handling different types of boundary conditions and averts solving complicated transcendental equations as demanded by eigenfunction expansion methods. Further works involving the Laplace transform have studied the permeable layered reaction diffusion problem in [41,42]. Solutions obtained in these works are restricted to two layers as well as obtaining the inverse Laplace transform numerically. In the same context, generalized integral transform techniques, for short GITT, are well-established hybrid approaches for solving diffusion and convection-diffusion problems, in which hybrid refers to the combination of classical analytical methods with modern computational tools aimed at accurate, robust, and low-cost solutions [43][44][45][46][47]. In the current work, we aim to extend, generalize, and merge results in [31,33,38,40,42] to solve certain nonhomogeneous diffusion problems in one-dimensional n-layered media. We use a new generalized integral transform recently introduced in [48]. The obtained solutions are applicable to more general linear nonhomogeneous diffusion equations, finite media consisting of arbitrary many layers, continuity and dispersive flow at the contact interfaces between sequal layers, and transitory boundary conditions of the arbitrary type at the inlet and outlet. To the best knowledge of the authors, analytical solutions verifying all the above mentioned conditions have not been previously reported in the literature which strongly motivates this current work.
In the remaining part of this introductory section, in Subsection 1.1, the multilayer diffusion problem is described, and then, it is reformulated as a sequence of one-layer diffusion problems having boundary conditions including given time-depending functions. Basic properties for theM ρ,m -transform that will be needed in this work are stated in Subsection 1.2. The remaining sections are constructed as follows: in Section 2, we discuss the solvability of a general linear nonhomogeneous one-layer diffusion problem with arbitrary time-varying data, using the M ρ,m -transform. Section 3 is devoted to our main multilayer diffusion problem, where in Subsection 3.1, we solve a two-layer problem to shed light on the basic idea by considering this simple case. Further, in Subsection 3.2, we return to benefit from the results obtained in Section 2 and Subsection 3.1 to solve the main multilayer diffusion problems (2)-(8) (see Subsection 1.1 below).

Mathematical Modeling for Nonhomogeneous n-Layer
Diffusion Systems. A one-dimensional diffusion problem in an n-layered permeable medium is set out as follows. Let be a finite partition of the interval ½α, β. In each subinterval ½x j−1 , x j , with j = 1, 2, ⋯, n, the component function φ j ðx, tÞ satisfies the partial differential equation (PDE) where d j ≥ 0, for all 1 ≤ j ≤ n, are the diffusion coefficients and with m ∈ ℤ + = f1, 2, 3,⋯g, ρ ∈ ℂ, ReðρÞ > 0. Here, the function-term λðt, τÞr j ðx, tÞ physically means the external source term that could be applied to the diffusion equation with r j ðx, tÞ depends on time and space while the other factor of the source term, i.e., λ, depends only on time. This last factor could be, for instance, a periodic-time magnetic source. The initial conditions (ICs) are assumed as The boundary conditions (BCs) are posited as (i) The outer BCs (at the inlet x = α and the outlet x = β) are general Robin boundary conditions as for all t ≥ 0, with ı, ι, ℓ, and l are constants satisfying |ı | + | ι | >0, |ℓ | +|l| > 0.
For appropriate given functions η 1 , ⋯, η n , ζ, and ξ, we are going to find an analytic solution of the problems (2)-(8) using the M ρ,m -generalized integral transform, introduced recently in [48]. Problems (2)-(8) can be reduced into the following sequence of one-layer diffusion problems.

2
Journal of Function Spaces (i) In the inlet layer, i.e., x ∈ ½x 0 , x 1 , (ii) In the interior layers, i.e., x ∈ ½x j−1 , x j , 2 ≤ j ≤ n − 1, (iii) In the outlet layer, i.e., x ∈ ½x n−1 , x n , Remark 1. Each of the initial boundary value problems (9)-(11) is a case of the one-layer nonhomogeneous diffusion problem that will be discussed in Section 2 below. Now, in view of the inner boundary conditions (7) and (8), the time-varying functions ζ j and ξ j for all 2 ≤ j ≤ n are subject to so that While the outer boundary data ζ 1 ðtÞ = ζðtÞ and ξ n ðtÞ = ξðtÞ are given in (5) and (6), respectively, the functions ζ j ð 2 ≤ j ≤ nÞ can be determined once we specify the functions ξ j ð1 ≤ j ≤ n − 1Þ. Hence, we have to find ξ j , 1 ≤ j ≤ n − 1. To do so, we should use the first matching condition (7).

Srivastava-Luo-Raina Generalized Integral Transform.
In [48], Srivastava et al. introduced the following generalized integral transform: for a continuous (or piecewise continuous) function φ on ½ 0, ∞Þ, where ρ ∈ ℂ ; Re ðρÞ ≥ 0 ; m ∈ ℤ + , s > 0 is the transform variable and τ > 0 is a parameter. The basic properties of the M ρ,m -transform are given in [48]. Next, we recall some of these properties, which are needed in the present work. Indeed, as introduced in [48] the M ρ,m -transform is closely related with the well-known integral transforms, the Laplace, natural, and Sumudu transforms. The Laplace transform is defined by So, from (14) and (15), we have the following duality relations: Setting ρ = 0 in (14), we recover the natural transform defined as (see [49,50])

Journal of Function Spaces
Thus, we have the following M ρ,m -ℕ-transform duality The Sumudu transform is defined by [51][52][53] Thus, Based on these dualities of the M ρ,m -transform (14) and these well-known integral transforms, it seems to be interesting to apply the M ρ,m -transform (14) in solving a variety of boundary and initial-boundary value problems. In this context, we recall the following results [48]: (i) Let φ ðnÞ ðtÞ be the n th -order t-derivative of the function φðtÞ and jφðtÞj ≤ Ke t/γ with K > 0, γ > 0. Then, where ℕ½φðtÞðs, τÞ is defined by (18). Using the duality (21) in (25), we find (ii) Again, using the dualities stated before a convolution formula for the M ρ,m -transform (14) can be obtained as follows. Here, the convolution for the Laplace transform will be considered; that is, for the functions φ and ψ, the convolution formula is given as If Φðs, τÞ = M ρ,m ½φðtÞðs, τÞ and Ψðs, τÞ = M ρ,m ½ψðtÞðs, τÞ, then-τΦðs, τÞΨðs, τÞ = τ Setting t 1 + t 2 = t in the last equality, one gets Here, changing of the integral order is used. Thus, using the duality of the M ρ,m and ℕ transforms (see (20)), we find Remark 2. If we put ρ = 0 in (30), the case being interesting later in our work, then we get (iii) Once again, ...again, using the dualities stated before an inversion formula of the M ρ,m -transform (14) is as long as the integral converges absolutely. In case, when ρ = 0 one obtains the following inversion formula of the natural transform (see Theorem 5.3 of [49]) The residue theorem (see, e.g., [54]) is usually used to calculate the contour integrals in (32) and (33).
Applying the M ρ,m -transform defined by (14), to (34), yields Using the duality of the M ρ,m -transform and the natural transform given by (21) and (25), Equation (38) can be reduced to where ℕ½φðtÞðs, τÞ is defined by (18).
then, (39) can be expressed as where Applying the variation of the parameter method to the nonhomogeneous equation (41) gives the general solution as where A and B are arbitrary invariants which can depend on s and τ. Differentiating (43) Transforming the boundary conditions (36) and (37), implies For simplicity, we set the following vector notations: Obviously, we have Substituting (43) and (44) into (45) and using the vector notation, we give the algebraic linear system where h·i is the usual dot product in ℝ 2 , and with Hðs, τÞ = ℕ½λðt, τÞξðtÞ and ζðtÞ and ξðtÞ are the boundary data given in (36) and (37), respectively. The  (50) is is the determinant of the coefficient matrix of system (50). Substituting the constants A and B into (43) gives which can be rewritten as where For further computation, we rewrite b θðx ; s, τÞ as Lemma 3. Let τ, s, x, y ∈ ℝ and L ∈ ℝ 2 . Then, Consequently, for each zero s of the function ΔðsÞ sinh ffiffiffiffiffiffiffiffiffiffiffiffi ðs/τdÞ p ðx − yÞ, one has where ψ is given by (57).
Proof. The first two conclusions of the lemma follow directly from the uniqueness theorem of the initial value problem for the second-order ordinary differential equations having constant coefficients.
It is easy to see that as functions in x, both sides of (61) are linear combinations of the functions cosh ffiffiffiffiffiffiffiffiffiffiffiffi ðs/τdÞ Journal of Function Spaces and sinh ffiffiffiffiffiffiffiffiffiffiffiffi ðs/τdÞ p x which are linearly independent solutions to the differential equation in (63). Thus, both sides of (61) solve this differential equation.
Moreover, in view of (57), both sides of (61) satisfy Hence, by the uniqueness theorem, (61) holds true. For each s * being a zero of the function ΔðsÞ sinh ffiffiffiffiffiffiffiffiffiffiffiffi ðs/τdÞ p ðx − yÞ, taking the limit in both sides of (61) as s ⟶ s * gives (62).
Applying Lemma 3, (56) and (59) respectively can be reduced to Next, in order to obtain the solution to the initial value problem (34)-(37), we apply the inversion formula (33) to (65) and (66). In doing so, we suppose that there are nonzero simple roots fs k g ∞ k=1 of ΔðsÞ. That is, Lemma 4. Suppose that (67) holds true. For each x, y ∈ ℝ and t, τ > 0, we get where with Proof. Let Applying the inversion formula (33), we find The last integral can be usually calculated by the residue theorem [54]. Hence, Recalling (67), each s k ðk = 1, 2,⋯Þ is a simple pole of e st/τ Rðs, τÞ. Therefore,

Journal of Function Spaces
We see that either as s tends to 0. Then, s = 0 is either a removable singular point or a simple pole of e st/τ Rðs, τÞ. Hence, substituting (74) and (75) in (73) gives the main conclusion of the Lemma, i.e., (69) and (70).☐ In view of (68) and (42), we have where Θ is defined by (69) and (70), and Hence, (66) can be rewritten as By the convolution formula (31), the inverse natural transform of (79) is That is, From Lemma 3, one has at s = 0 and s = s k ðk = 1, 2,⋯Þ the zeros of ΔðsÞ sinh ffiffiffiffiffiffiffiffi s/τd p ð x − yÞ. That results in The first conclusion is obvious when Θ 0 = 0 in (70). Thus, (81) can be simplified as where b θðx ; s, τÞ is given in (66). Now, we can obtain the solution φðx, tÞ of Problem (34)-(37) by operating the inversion formula (33) in (85). In doing so, we need the following lemma.

Illustrative Examples.
Here, we discuss two illustrative test cases to show the accuracy and effectiveness of our technique.
Example 1. Heat equation with zero temperatures at finite ends.

Multilayer Nonhomogeneous Diffusion System
Here, we are seeking the solution of our main problem defined in (2)-(8), which was converted into a sequence of initial boundary value problems (9)- (11). For the convenience of the reader and in order to draw the full picture in an easy way, we start with solving the bilayer diffusion problem in the following subsection; then, we move to the general case in Section 3.2.