Modeling diffusion and convection in a complex system

[cg78ithaca]

I am trying to come up with an analytical solution (even as a infinite series etc.) for the following diffusion-convection problem.

A thin layer of gel (assumed rectangular) is in direct contact with a liquid layer (perfusate) flowing with velocity v in the x direction (left to right) just underneath it. On the side of the gel not in contact with the perfusate a solid drug is sprayed. As the solid drug dissolves into the gel layer, it slowly diffuses across (in the z+ direction) into the liquid layer, from where it is transported by convection. The problem is finding an analytical expression for the exit concentration, that is, the concentration of drug in the perfusate layer at the last point on the x-axis where the gel and perfusate are in contact (x = L_x).

A diagram illustrating the system is shown below. Although the system is 3D, the solid drug is assumed to be sprayed uniformly on the (x,y,z) = (x,0,z) face of the gel phase parallelepiped so that the solution will be uniform in y at all times. This reduces the problem to a 2D one:

Several simplifications are then made. The dissolution of the solid drug into the gel is assumed to be instantaneous and by way of a superficial layer which is permanently saturated. That is, the region of the gel between z = 0 and z = δ_z is permanently at C_s (saturation concentration) and diffusion occurs from this layer into the depth of the gel.

Furthermore, since diffusion in the x-axis progresses much more slowly than convection, it is neglected in both the gel and perfusate phase. This reduces the problem to two differential equations, one for the gel phase, with their appropriate boundary conditions. C refers to the gel phase and K refers to the perfusate phase:
$$
\frac{\partial{}C_{gel}}{\partial{}t}=D_g\frac{{\partial{}}^2}{\partial{}z^2}C_{gel}
$$ and one for the perfusate phase:
$$
\frac{\partial{}K}{\partial{}t}=D_p\frac{{\partial{}}^2K}{\partial{}z^2}-v\frac{\partial{}K}{\partial{}x}
$$
Introducing the notation ζ= z - δ_z and l_z = L_zg - δ_z, the gel equation becomes:
$$
\frac{\partial{}C_{\zeta{}}}{\partial{}t}=D_g\frac{{\partial{}}^2}{\partial{}{\zeta{}}^2}C_{\zeta{}}
$$

Steady-State Solution

Attempting the method of separation of variables in the case of the gel equation for example leads generally to a sum of exponentials in t, and requiring that this does not diverge as t --> Inf means that they are all of the form exp(-kt), for k >0. That means that as t --> Inf the separation of variables solution will converge to zero, which means that the steady-state solution will have to be added to that (equivalent to adding the case k = 0). That is ok since the differential equations are linear and the steady-state is a solution to them. Thus the initial problem is finding the steady-state solutions, governed by the equations:
$$
0=\frac{\partial{}C_{\zeta{},\infty{}}}{\partial{}t}=D_g\frac{{\partial{}}^2}{\partial{}{\zeta{}}^2}C_{\zeta{},\infty{}}
$$ subject to the boundary condition at ζ = 0:
$$
C_{\zeta{},\infty{}}\left(x,\zeta{}=0\right)=C_S,\ \forall{}x
$$

For the perfusate phase (note that z restarts at zero at the gel/perfusate interface in K):
$$
0=\frac{\partial{}K}{\partial{}t}=D_p\frac{{\partial{}}^2K}{\partial{}z^2}-v\frac{\partial{}K}{\partial{}x}
$$ subject to the boundary conditions at x = 0 requiring no drug in the incoming perfusate and a closed boundary condition at z = L_zp (bottom of perfusate):
$$
K_{\infty{}}\left(x=0,z\right)=0,\forall{}z
$$ $$\frac{\partial{}K_{\infty{}}}{\partial{}z}\left(x,z=0\right)=0,\forall{}x
$$
At the interface the two equations are connected via the boundary conditions there, which require equipartition of drug between the layers.
$$
\alpha{}C_{\zeta{},\infty{}}(x,l_z)=K_{\infty{}}(x,0)=f(x),\forall{}x
$$ and conservation of mass at the interface, which requires that the fluxes out of the gel and into the perfusate are equal:
$$
D_g\frac{\partial{}C_{\zeta{},\infty{}}}{\partial{}\zeta{}}\left(x,l_z\right)=D_p\frac{\partial{}K_{\infty{}}}{\partial{}z_p}\left(x,0\right)=D_p g(x),\forall{}x
$$ Note I'm introducing two functions here, f(x) and g(x), referring to the concentration at the interface along x and the flux at the interface along x, respectively.

These conditions will also have to hold for all t, not just at steady-state.

The steady-state solution for the gel phase is simple to derive since it is a first-order polynomial, subject to the respective boundary conditions:
$$
C_{\zeta{},\infty{}}\left(x,\zeta{}\right)=C_S+\left\{\frac{1}{\alpha{}l_z}f\left(x\right)-\frac{C_S}{l_z}\right\}\zeta{}
$$ which imposes the following requirement on the flux across the interface:
$$
g\left(x\right)=\frac{D_g}{l_zD_p}\left\{\frac{1}{\alpha{}}f\left(x\right)-C_S\right\}
$$
The steady-state solution for the perfusate layer right under the gel is more difficult to arrive at, but the requirement of zero flux at z = L_zp and fixed slope and intercept at z = 0 via f(x) and g(x) lead to a special quantization requirement for the standing waves that make up the general solution (k = D_p/v):
$$
K_{\infty{}}\left(z,t\right)=\alpha{}C_S-\sum_{n=0,1,2…}^{\infty{}}{\epsilon{}}_n\left\{\frac{\sqrt{k}}{{\lambda{}}_n}\sin{\left(\frac{{\lambda{}}_n}{\sqrt{k}}z\right)}+\alpha{}\frac{l_zD_p}{D_g}\cos{\left(\frac{{\lambda{}}_n}{\sqrt{k}}z\right)}\right\}\exp{\left(-{{\lambda{}}_n}^2t\right)}
$$ with
$$
{\lambda{}}_n=\frac{D_g}{l_zD_p}\frac{\sqrt{k}}{\alpha{}}\cot{\left(\frac{{\lambda{}}_n}{\sqrt{k}}L_{zp}\right)},\
\frac{{\lambda{}}_n}{\sqrt{k}}L_{zp}\in{}\left(n\pi{},\left(n+\frac{1}{2}\right)\pi{}\right),n=0,1,2…
$$ $${\epsilon{}}_n=\frac{2\alpha{}C_S}{{\left(\alpha{}\frac{l_zD_p}{D_g}\right)}^2\frac{{{\lambda{}}_n}^2}{k}L_{zp}+\alpha{}\frac{l_zD_p}{D_g}+L_{zp}}
$$
Since cot() is a surjective function on the real axis for each of its periods, there will be a unique solution for each period of cot() of the lambda-equation above, as shown below in a plot of λ_nL_zp/√k – nπ in units of π vs. n for some typical values of the parameters:

These standing waves turn out to form an orthogonal set, which allows for determining the amplitudes much as one would for a Fourier series, and a plot of the ε_n is shown below:

They converge nicely as n increases, which allows for determining the general solution K(x,z,t->∞) to arbitrary precision, e.g. as a function of z at some length x > 0:

This has the expected features, an oblique asymptote at z = 0 exactly matching the g(x) requirement and and horizontal asymptote at z = L_zp.

A plot of f(x)/α to show the concentration profile at the interface is below (i.e. K(x,z=0,t->∞)/α as a function of x:

Finally, while for x < 0 the concentration profile is clearly K(x,z,t) = 0 (since there is no drug in the incoming perfusate), past x > L_x it is a simple case of diffusion in the z axis subject to closed boundary conditions at z = 0 and z = L_zp and convection in the x direction, which solves to (by Fourier series as a sum of cosines given the boundary conditions):

$$
K_{\infty{}}\left(x,z\right)=\frac{1}{2}B_0+\sum_{m=1,2,..}^{\infty{}}B_m\cos{\left(\frac{m\pi{}}{L_{zp}}z\right)}\exp{\left(-k{\left(\frac{m\pi{}}{L_{zp}}\right)}^2\left(x-L_x\right)\right)},x>L_x
$$ with
$$
B_0=2\alpha{}C_S-\frac{2k}{L_{zp}}\sum_{n=0,1,2…}^{\infty{}}\frac{{\epsilon{}}_n}{{{\lambda{}}_n}^2}\exp{\left(-{{\lambda{}}_n}^2L_x\right)}
$$ $$
B_m=2L_{zp}\sum_{n=0,1,2…}^{\infty{}}\frac{{\epsilon{}}_n}{{\pi{}}^2m^2-\frac{{{\lambda{}}_n}^2}{k}{L_{zp}}^2}\exp{\left(-{{\lambda{}}_n}^2L_x\right)}
$$

Time-Dependent Solution

With the steady-state solutions available, it is possible now to express the general time-dependent solution in the gel phase as a component that converges to zero as t->∞ plus the steady-state:
$$
C_{\zeta{}}\left(x,\zeta{},t\right)=C_S+\left\{\frac{1}{\alpha{}l_z}f\left(x\right)-\frac{C_S}{l_z}\right\}\zeta{}-u\left(x,\zeta{},t\right)
$$
The u(x,ζ,t) function is the solution to:
$$
\frac{\partial{}u}{\partial{}t}=D_g\frac{{\partial{}}^2u}{\partial{}{\zeta{}}^2}
$$ with the initial condition:
$$
u\left(x,\zeta{},0\right)=\left\{\begin{array}{l}0\ if\
\zeta{}=0,x\in{}\left[0,L_x\right] \\
C_S-\zeta{}\frac{l_zD_p}{D_g}\sum_{n=0,1,2…}^{\infty{}}{\epsilon{}}_n\exp{\left(-{{\lambda{}}_n}^2x\right)}\
for\ \zeta{}\in{}\left(0,l_z\right),x\in{}\left[0,L_x\right]\end{array}\right.
$$ boundary condition at ζ = 0:
$$
u\left(x,0,t\right)=0,\ \forall{}x\in{}\left[0,L_x\right],t\geq{}0
$$ and interface matching conditions:
$$
u\left(x,l_z,t\right)=C_S-\frac{1}{\alpha{}}f\left(x,t\right)-\frac{l_zD_p}{D_g}\sum_{n=0,1,2…}^{\infty{}}{\epsilon{}}_n\exp{\left(-{{\lambda{}}_n}^2x\right)},\forall{}x\in{}\left[0,L_x\right],t\geq{}0
$$ $$
\frac{\partial{}}{\partial{}\zeta{}}u\left(x,l_{z},t\right)=-\frac{D_p}{D_g}\left\{g\left(x,t\right)+\sum_{n=0,1,2…}^{\infty{}}{\epsilon{}}_n\exp{\left(-{{\lambda{}}_n}^2x\right)}\right\},\forall{}x\in{}\left[0,L_x\right],t\geq{}0
$$
Note that now f and g have become functions of x and t since we are interested in the general time-domain solution, not just the steady-state.

This problem can also be solved by separation of variables, with the solution of the form:
$$
u\left(x;\zeta{},t\right)=\int_{-\infty{}}^{\infty{}}\Lambda{}(x,\kappa{})\sin{\left(\frac{\kappa{}}{\sqrt{D_g}}\zeta{}\right)}\exp{\left(-{\kappa{}}^2t\right)}d\kappa{}
$$ where Λ is a function determined by:
$$
\int_{-\infty{}}^{\infty{}}\Lambda{}(x,\kappa{})\sin{\left(\frac{\kappa{}}{\sqrt{D_g}}l_z\right)}d\kappa{}=u\left(x;\zeta{},0\right)=\left\{\begin{array}{l}0\
if\ \zeta{}=0 \\
C_S-\zeta{}\frac{l_zD_p}{D_g}\sum_{n=0,1,2…}^{\infty{}}{\epsilon{}}_n\exp{\left(-{{\lambda{}}_n}^2x\right)}\
otherwise\end{array}\right.
$$
Therefore f(x,t) and g(x,t) are:
$$
f\left(x;t\right)=\alpha{}C_S-\alpha{}l_z\frac{D_p}{D_g}\sum_{n=0,1,2…}^{\infty{}}{\epsilon{}}_n\exp{\left(-{{\lambda{}}_n}^2x\right)}-\alpha{}\int_{-\infty{}}^{\infty{}}\Lambda{}(x,\kappa{})\sin{\left(\frac{\kappa{}}{\sqrt{D_g}}l_z\right)}\exp{\left(-{\kappa{}}^2t\right)}d\kappa{}
$$ $$
g\left(x,t\right)=-\sum_{n=0,1,2…}^{\infty{}}{\epsilon{}}_n\exp{\left(-{{\lambda{}}_n}^2x\right)}-\frac{\sqrt{D_g}}{D_p}\int_{-\infty{}}^{\infty{}}\Lambda{}(x,\kappa{})\kappa{}\cos{\left(\frac{\kappa{}}{\sqrt{D_g}}l_z\right)}\exp{\left(-{\kappa{}}^2t\right)}d\kappa{}
$$
In the spirit of the steady-state solution, it is possible that κ only takes on discrete values, in which case either the Λ function can be seen as a sum of delta functions, or alternatively as a sum:
$$
u\left(x;\zeta{},t\right)=\sum_j{\Lambda{}}_j(x)\sin{\left(\frac{{\kappa{}}_j}{\sqrt{D_g}}\zeta{}\right)}\exp{\left(-{\kappa{}}_j^2t\right)}
$$ with Λ subject to:
$$
\sum_j{\Lambda{}}_j(x)\sin{\left(\frac{{\kappa{}}_j}{\sqrt{D_g}}\zeta{}\right)}=u\left(x;\zeta{},0\right)=\left\{\begin{array}{l}0\
if\ \zeta{}=0 \\
C_S-\zeta{}\frac{l_zD_p}{D_g}\sum_{n=0,1,2…}^{\infty{}}{\epsilon{}}_n\exp{\left(-{{\lambda{}}_n}^2x\right)}\
otherwise\end{array}\right.
$$
Therefore f(x,t) and g(x,t) are:
$$
f\left(x;t\right)=\alpha{}C_S-\alpha{}l_z\frac{D_p}{D_g}\sum_{n=0,1,2…}^{\infty{}}{\epsilon{}}_n\exp{\left(-{{\lambda{}}_n}^2x\right)}-\alpha{}\sum_j{\Lambda{}}_j(x)\sin{\left(\frac{{\kappa{}}_j}{\sqrt{D_g}}\zeta{}\right)}\exp{\left(-{\kappa{}}_j^2t\right)}
$$ $$
g\left(x,t\right)=-\sum_{n=0,1,2…}^{\infty{}}{\epsilon{}}_n\exp{\left(-{{\lambda{}}_n}^2x\right)}-\frac{\sqrt{D_g}}{D_p}\sum_j{\Lambda{}}_j(x){\kappa{}}_j\cos{\left(\frac{{\kappa{}}_j}{\sqrt{D_g}}\zeta{}\right)}\exp{\left(-{\kappa{}}_j^2t\right)}
$$
As before, the restrictions on κ would arise from matching the two equations at the interface.

 

[cg78ithaca]

continued below:

Similarly for the time-dependent perfusate equation, writing the time-dependent solution as a component which converges to zero as t->∞ plus the steady-state:
$$
K\left(x,z,t\right)=u\left(x,z,t\right)+K_{\infty{}}\left(x,z\right)
$$ The function u satisfies a PDE of the type:
$$
\frac{\partial{}u}{\partial{}t}=D_p\frac{{\partial{}}^2u}{\partial{}z^2}-v\frac{\partial{}u}{\partial{}x}
$$ which can be simplified by transforming to a frame of reference moving with velocity v, i.e. introducing ξ(x,t) = x – vt, w(ξ(x,t),z,t) = u(x,z,t):
$$
\frac{\partial{}w}{\partial{}\xi{}}=D_p\frac{{\partial{}}^2w}{\partial{}z^2}
$$ which is now subject to the initial condition (note that unlike the gel case, in the perfusate 0 < x < ∞):
$$
w\left(\xi{},z,0\right)=\left\{\begin{array}{l}0\ if\ \xi{}\leq{}0 \\
-\alpha{}C_S+\sum_{n=0,1,2…}^{\infty{}}{\epsilon{}}_n\left\{\frac{\sqrt{k}}{{\lambda{}}_n}\sin{\left(\frac{{\lambda{}}_n}{\sqrt{k}}z\right)}+\alpha{}\frac{l_zD_p}{D_g}\cos{\left(\frac{{\lambda{}}_n}{\sqrt{k}}z\right)}\right\}\exp{\left(-{{\lambda{}}_n}^2\xi{}\right)}if\
0<\xi{}\leq{}L_x \\
-\frac{1}{2}B_0-\sum_{m=1,2,..}^{\infty{}}B_m\cos{\left(\frac{m\pi{}}{L_{zp}}z\right)}\exp{\left(-k{\left(\frac{m\pi{}}{L_{zp}}\right)}^2\left(\xi{}-L_x\right)\right)}if\
\xi{}>L_x\end{array}\right.
$$ The boundary conditions are:
$$
w\left(-vt,z,t\right)=u\left(0,z,t\right)=0,\forall{}z\in{}\left[0,L_{zp}\right],\forall{}t>0
$$ $$
w_z\left(\xi{},L_{zp},t\right)=u_z\left(x,L_{zp},t\right)=0,\forall{}\xi{}\in{}R,\forall{}t>0
$$ with the interface boundary condition equations:
$$
w\left(\xi{},0,t\right)=\left\{\begin{array}{l}0\ if\ \xi{}<-vt,\forall{}t>0 \\
f\left(\xi{}+vt,t\right)-\alpha{}C_S+\alpha{}\frac{l_zD_p}{D_g}\sum_{n=0,1,2…}^{\infty{}}{\epsilon{}}_n\exp{\left(-{{\lambda{}}_n}^2(\xi{}+vt)\right)}
\\
not\ restricted,\forall{}\
\xi{}>L_x-vt,\forall{}t>0\end{array}\right.if-vt\leq{}\xi{}\leq{}L_x-vt,\forall{}t>0
$$ $$
w_z(\xi{},0,t)=\left\{\begin{array}{l}0\ if\
\xi{}<-vt,\forall{}t>0 \\
g\left(\xi{}+vt,t\right)+\sum_{n=0,1,2…}^{\infty{}}{\epsilon{}}_n\exp{\left(-{{\lambda{}}_n}^2(\xi{}+vt)\right)},-vt\leq{}\xi{}\leq{}L_x-vt
\\
0,\forall{}\xi{}>L_x-vt,\forall{}t>0\end{array},\forall{}t>0\right.
$$
Since t does not appear at all in the PDE:
$$
\frac{\partial{}w}{\partial{}\xi{}}=D_p\frac{{\partial{}}^2w}{\partial{}z^2}
$$ the time component of w(ξ,z,t) is completely separable, and generally
spanning the space of all well-behaved functions T_q(t):
$$
w\left(\xi{},z,t\right)=\sum_q{\psi{}}_q\left(z,\xi{}\right)T_q(t)
$$ Then:
$$
\sum_q{\psi{}}_q\left(z,\xi{}\right)T_n(0)=\left\{\begin{array}{l}0\ if\
\xi{}\leq{}0 \\
-\alpha{}C_S+\sum_{n=0,1,2…}^{\infty{}}{\epsilon{}}_n\left\{\frac{\sqrt{k}}{{\lambda{}}_n}\sin{\left(\frac{{\lambda{}}_n}{\sqrt{k}}z\right)}+\alpha{}\frac{l_zD_p}{D_g}\cos{\left(\frac{{\lambda{}}_n}{\sqrt{k}}z\right)}\right\}\exp{\left(-{{\lambda{}}_n}^2\xi{}\right)}if\
0<\xi{}\leq{}L_x \\
-\frac{1}{2}B_0-\sum_{m=1,2,..}^{\infty{}}B_m\cos{\left(\frac{m\pi{}}{L_{zp}}z\right)}\exp{\left(-k{\left(\frac{m\pi{}}{L_{zp}}\right)}^2\left(\xi{}-L_x\right)\right)}if\
\xi{}>L_x\end{array}\right.
$$ which is an expansion of the initial condition in the basis set ψ_q.
The boundary conditions are:
$$
\sum_qT_q\left(t\right){\psi{}}_q\left(-vt,z\right)=w\left(-vt,z,t\right)=0,\forall{}z\in{}\left[0,L_{zp}\right],\forall{}t>0
$$ $$\sum_qT_q\left(t\right){\psi{}}_{q,z}\left(\xi{},L_{zp}\right)=w_z\left(\xi{},L_{zp},t\right)=0,\forall{}\xi{}\in{}R,\forall{}t>0
$$
with the connection equations:
$$
\sum_qT_q\left(t\right){\psi{}}_q\left(\xi{},0\right)=\left\{\begin{array}{l}0\
if\ \xi{}<-vt,\forall{}t>0 \\
f\left(\xi{}+vt,t\right)-\alpha{}C_S+\alpha{}\frac{l_zD_p}{D_g}\sum_{n=0,1,2…}^{\infty{}}{\epsilon{}}_n\exp{\left(-{{\lambda{}}_n}^2(\xi{}+vt)\right)}
\\
not\ restricted,\forall{}\
\xi{}>L_x-vt,\forall{}t>0\end{array}\right.if-vt\leq{}\xi{}\leq{}L_x-vt,\forall{}t>0
$$ $$
\sum_qT_q(t){\psi{}}_{q,z}(\xi{},0)=\left\{\begin{array}{l}0\ if\
\xi{}<-vt,\forall{}t>0 \\
g\left(\xi{}+vt,t\right)+\sum_{n=0,1,2…}^{\infty{}}{\epsilon{}}_n\exp{\left(-{{\lambda{}}_n}^2(\xi{}+vt)\right)},-vt\leq{}\xi{}\leq{}L_x-vt
\\
0,\forall{}\xi{}>L_x-vt,\forall{}t>0\end{array},\forall{}t>0\right.
$$ To figure out what this ψ_q basis set might be, take the PDE:
$$
\frac{\partial{}{\psi{}}_q}{\partial{}\xi{}}=D_p\frac{{\partial{}}^2{\psi{}}_q}{\partial{}z^2}
$$ which solves generally by separation of variables as seen before to:
$$
{\psi{}}_q\left(z,\xi{}\right)=\int_0^{\infty{}}\left\{A_q(\rho{})\sin{\left(\frac{\rho{}}{\sqrt{D_p}}z\right)}+B_q(\rho{})\cos{\left(\frac{\rho{}}{\sqrt{D_p}}z\right)}\right\}\exp{\left(-{\rho{}}^2\xi{}\right)}d\rho{}
$$
Note that ψ_q(z,ξ) -> 0 as ξ -> ∞ and that a constant function also solves the equation above; let that constant function solution be ψ₀.
Expressing ψ_q as the integral shown above may actually be overkill, as I may simply have a discrete set of ρ's and each ψ_n simply the function corresponding to each ρ_n, or some linear combination thereof.
Then:
$$
w\left(\xi{},z,t\right)={\psi{}}_0+\sum_{q\not=0}{\psi{}}_q\left(z,\xi{}\right)T_q(t)
$$ The initial condition requires, for ξ >= 0:
$$
{\psi{}}_0+\int_0^{\infty{}}\exp{\left(-{\rho{}}^2\xi{}\right)}\sum_{q\not=0}T_q(0)\left\{A_q(\rho{})\sin{\left(\frac{\rho{}}{\sqrt{D_p}}z\right)}+B_q(\rho{})\cos{\left(\frac{\rho{}}{\sqrt{D_p}}z\right)}\right\}d\rho{}=\left\{\begin{array}{l}-\alpha{}C_S+\sum_{n=0,1,2…}^{\infty{}}{\epsilon{}}_n\left\{\frac{\sqrt{k}}{{\lambda{}}_n}\sin{\left(\frac{{\lambda{}}_n}{\sqrt{k}}z\right)}+\alpha{}\frac{l_zD_p}{D_g}\cos{\left(\frac{{\lambda{}}_n}{\sqrt{k}}z\right)}\right\}\exp{\left(-{{\lambda{}}_n}^2\xi{}\right)}if\
0<\xi{}\leq{}L_x \\
-\frac{1}{2}B_0-\sum_{m=1,2,..}^{\infty{}}B_m\cos{\left(\frac{m\pi{}}{L_{zp}}z\right)}\exp{\left(kL_x{\left(\frac{m\pi{}}{L_{zp}}\right)}^2\right)}\exp{\left(-k{\left(\frac{m\pi{}}{L_{zp}}\right)}^2\xi{}\right)}if\
\xi{}>L_x\end{array}\right.
$$ Taking now the limit of ξ -> ∞, it follows that:
$$
{\psi{}}_0=-\frac{1}{2}B_0
$$ Now substituting τ = ρ² > 0:
$$
\int_0^{\infty{}}\exp{\left(-\tau{}\xi{}\right)}\sum_q\frac{T_q(0)}{2\sqrt{\tau{}}}\left\{A_q^*(\tau{})\sin{\left(\sqrt{\frac{\tau{}}{D_p}}z\right)}+B_q^*(\tau{})\cos{\left(\sqrt{\frac{\tau{}}{D_p}}z\right)}\right\}d\tau{}=\left\{\begin{array}{l}\frac{1}{2}B_0-\alpha{}C_S+\sum_{n=0,1,2…}^{\infty{}}{\epsilon{}}_n\left\{\frac{\sqrt{k}}{{\lambda{}}_n}\sin{\left(\frac{{\lambda{}}_n}{\sqrt{k}}z\right)}+\alpha{}\frac{l_zD_p}{D_g}\cos{\left(\frac{{\lambda{}}_n}{\sqrt{k}}z\right)}\right\}\exp{\left(-{{\lambda{}}_n}^2\xi{}\right)}if\
0<\xi{}\leq{}L_x \\
-\sum_{m=1,2,..}^{\infty{}}B_m\cos{\left(\frac{m\pi{}}{L_{zp}}z\right)}\exp{\left(kL_x{\left(\frac{m\pi{}}{L_{zp}}\right)}^2\right)}\exp{\left(-k{\left(\frac{m\pi{}}{L_{zp}}\right)}^2\xi{}\right)}if\
\xi{}>L_x\end{array}\right.
$$ Now this means that the RHS is the Laplace transform of the sum over q, and
therefore that:
$$
\sum_q\frac{T_q(0)}{2\sqrt{\tau{}}}\left\{A_q^*(\tau{})\sin{\left(\sqrt{\frac{\tau{}}{D_p}}z\right)}+B_q^*(\tau{})\cos{\left(\sqrt{\frac{\tau{}}{D_p}}z\right)}\right\}=L^{-1}\{RHS\}(\tau{})
$$
This brings me to the problem I posted here at
https://www.physicsforums.com/threa...sform-of-a-piecewise-defined-function.904387/
which revolves around determining the inverse Laplace transform of a piecewise-defined function as the one above is.