The diffusion equation with time-dependent boundary condition

[Corribus]

TL;DR Summary

Trying to use Laplace Transform method to solve a 1D diffusion problem with a time-dependent boundary condition. Looking for help identifying my error.

Hi everyone,

I am trying to solve the 1 dimensional diffusion equation over an interval of 0 < x < L subject to the boundary conditions that C = kt at x = 0 and C = 0 at x = L. k is a constant. The diffusion equation is

[tex]\frac{dC}{dt}=D\frac{d^2C}{dx^2}[/tex]

I am using the Laplace transform method, which I have applied successfully to solve similar diffusion problems.

Briefly, when applying the Laplace transform and solving the subsidiary equation with the stated boundary conditions, I arrive at the solution:

[tex]\bar C=\frac{k \sinh{q(L-x)}}{p^2 \sinh{qL}}[/tex]

Where q^2 = p/D and p is the transform variable.

To get the inverse transform, I’m using the partial fraction approach, basically expanding the hyperbolic sine functions as infinite products, finding the roots, and going from there. A method to do this is laid out in Crank’s book Mathematics of Diffusion and the method almost always works for me. But not this time.

Anyway, the solution I’m getting using this approach is:

[tex]C=kt \frac{L-x}{L}+k\frac{(L-x)^3}{8LD}+\frac{2kL^2}{D\pi^3}\sum_{n=1}^{\infty}\frac{1}{n^3}\sin{\left( \frac{n \pi x}{L} \right)}\exp{\left( -\frac{Dn^2 \pi^2 t}{L^2}\right)}[/tex]

Unfortunately the second term of this solution does not satisfy the boundary conditions, so I’ve done something wrong.

I did solve the problem a different way using a heat transfer book as a guide and the correct solution seems to be:

[tex]C=kt \frac{L-x}{L}-\frac{kL^2}{D}\left(\frac{x}{3L}-\frac{x^2}{2L^2}+\frac{x^3}{6L^3}\right)+\frac{2kL^2}{D\pi^3}\sum_{n=1}^{\infty}\frac{1}{n^3}\sin{\left( \frac{n \pi x}{L} \right)}\exp{\left( -\frac{Dn^2 \pi^2 t}{L^2}\right)}[/tex]

The good news is I can get a solution to my diffusion problem, but I’d still like to know what I’m doing wrong with my inverse transform, because this is usually my method of choice. It has something to do with the triple zero root in the denominator of the subsidiary eqn solution because that gives rise to the second term that is wrong. But I’ve gone over my work a zillion times and I can’t find my error.

If anyone has any suggestions, I’d really appreciate it. Happy to supply my intermediate math if that’s helpful for troubleshooting.

 

[anuttarasammyak]

Hi. I am curious to know how you get sinh formula by Laplace transformation.

 

[Corribus]

Hi. I am curious to know how you get sinh formula by Laplace transformation.

Thanks for your interest. First, I realize I forgot to state that the initial condition is C = 0 for the region 0 < x < L.

The solution to the subsidiary equation after taking the Laplace transform and applying the initial condition is:

[tex]\bar C = A \exp{\left( qx \right)} + B \exp{\left( -qx \right)}[/tex]

Transforming the boundary conditions leads to:

[tex] \bar C = \frac{k}{p^2}, x = 0[/tex]
[tex] \bar C = 0, x = L[/tex]

Using these conditions I solve for A and B, which have exponential terms that can be combined to form sinh functions. Actually the appearance of sinh and cosh functions is common in Laplace transformations of diffusion equations involving 1D diffusion through membranes with different surface boundary conditions.

 

[pasmith]

Using the hyperbolic sum formulae one can always rearrange [itex]u(x) = Ae^{kx} + Be^{-kx}[/itex] as [tex]C\cosh k(x - x_0) + D \sinh k(x - x_0).[/tex] One form has two arbitrary constants and the other has three, so you can choose one of them to your advantage. In initial value problems this enables you to just write down the solution from the initial values ([itex]C = u(x_0)[/itex], [itex]D = \frac1k u'(x_0)[/itex]) and in boundary value problems you take [itex]x_0[/itex] to be the boundary where the solution vanishes so that [itex]C = 0[/itex] and you can easily find [tex]
D = \frac{u(x_1)}{\sinh k(x_1-x_0)}.[/tex] The hyperbolic functions also have the advantage that one of them is even and the other is odd. Indeed only on an unbounded interval are the exponentials preferable, since finiteness at [itex]+\infty[/itex] requires [itex]A = 0[/itex] and finiteness at [itex]-\infty[/itex] requires [itex]B = 0[/itex].

The good news is I can get a solution to my diffusion problem, but I’d still like to know what I’m doing wrong with my inverse transform, because this is usually my method of choice. It has something to do with the triple zero root in the denominator of the subsidiary eqn solution because that gives rise to the second term that is wrong. But I’ve gone over my work a zillion times and I can’t find my error.

If anyone has any suggestions, I’d really appreciate it. Happy to supply my intermediate math if that’s helpful for troubleshooting.

We can't really help you if you don't show your working.

The easiest method probably is to calculate the residue at the origin as [tex]
\frac12 \left.\frac{d^2}{dp^2} (p^3 \bar C(p) )\right|_{p = 0}.[/tex]

 

[anuttarasammyak]

Summary:: Trying to use Laplace Transform method to solve a 1D diffusion problem with a time-dependent boundary condition. Looking for help identifying my error.
I am trying to solve the 1 dimensional diffusion equation over an interval of 0 < x < L subject to the boundary conditions that C = kt at x = 0 and C = 0 at x = L. k is a constant.

Matter keep pouring out from the Origin and diffuse but C=0 at x=L always. Where the matter goes ?

In order to explain it, I think there is a twin source at x=2L but minus sign C=-kt. You see the right answer shows odd property for x=L and a minus source at X=2L.

The sencon term of RHS is written as
[tex]-\frac{k}{6LD}x(x-L)(x-2L)[/tex]
It is obvious for the first and the third terms.

 

[Corribus]

@pasmith

That's really interesting about expressing the subsidiary solution in terms of cosh and sinh. I hadn't thought of that, but I can see how that saves a lot of time. Thanks for pointing it out.

The residue method that you mention is basically what I did, so it's good to have the confirmation that this is the right strategy. Clearly I'm doing something wrong with my derivatives, though. I will transcribe my working into LaTex when I get a chance, hopefully in the next day or two.

@anuttaraysamyak

The condition C = 0 at x = L would be used if the volume on the other side of the membrane is large (infinite) and well-stirred, so that any mass that transfers through the film becomes infinitely dilute instantaneously. It's a decent approximation but not really realistic as the time gets large or as the volume on the other side of the membrane is small. You can still calculate the amount of mass that passes through the interface at x = L at any time by taking the derivative of the mass flux at the interface and then integrating over time. The boundary condition can be adjusted to reflect a well-stirred finite volume, which is a little more realistic. But this is more of an intellectual exercise for me to better understand solution methods than for any practical purpose.

 

[Corribus]

Alright, I took me a while to transcribe all the equations to Tex, but below is my working of this problem. I tried to make sure I didn’t introduce any errors into the transcription. I have provided my whole solution, not just the parts I know to be problematic. I have indicated in bold text below where my work on the term that is causing me difficulty begins, in the event it is helpful to start there.

My work:

Starting with by subsidiary solution:

[tex] \bar C=\frac{k \sinh{q(L-x)}}{p^2 \sinh{qL}}[/tex]

With q2 = p/D. For shorthand I define the numerator as f(p) and the denominator as g(p).

First, I find the roots of g(p). The roots of p^2 are 0 and 0. The roots an of sinh(qL) in p are

[tex] a_n =\frac {i \sqrt{D} n \pi}{L}, n = 0, 1, 2 \ldots [/tex]

So, there is actually a triple root at 0.

For each nonzero (single) root an there is a term in Cx,t

[tex]\frac {f(a_n)}{g\prime(a_n)} e^{a_n t}[/tex]

Which after a lot of steps gives rise to the summation term in the solution provided in my earlier post, which I believe to be correct.

For each multiple root in g(p) of order m and value b, there is a term in the solution Cx,t according to

[tex]\sum_{s=0}^{m-1} \left[ \frac{d^s}{dp^s} \left\{ \frac {(p-b)^m f(p)}{g(p)} \right\} \right]_{p=b} \frac{t^{m-s-1}}{s!(m-s-1)!}e^{bt}[/tex]

In the case here, m = 3 and b = 0, so substituting in f(p) and g(p) above, there is a term in the solution Cx,t:

[tex]k \sum_{s=0}^{2} \left[ \frac{d^s}{dp^s} \left\{ \frac {p \sinh {q(L-x)}}{\sinh{qL}} \right\} \right]_{p=0} \frac{t^{3-s-1}}{s!(3-s-1)!}e^{bt}[/tex]

I take the different values of s in order. First, s = 0. Based on the previous equation, there is a term:

[tex] \left[ \frac {p \sinh {q(L-x)}}{\sinh{qL}} \right]_{p=0} \frac{t^{2}}{2}[/tex]

Which obviously equals zero.

Next, s = 1. There is a term:

[tex]kt \left[ \frac{d }{dp} \left\{ \frac {p \sinh {q(L-x)}}{\sinh{qL}} \right\} \right]_{p=0}[/tex]

I use the product rule for the derivative:

[tex]kt \left[ p \frac{d }{dp} \left\{ \frac {\sinh {q(L-x)}}{\sinh{qL}} \right\} + \left\{ \frac {\sinh {q(L-x)}}{\sinh{qL}} \right\} \right]_{p=0}[/tex]

From here forward, for shorthand, I let A = (L-x)/√D and B = L/√D, so this term rearranges to:

[tex]kt \left[ p \frac{d }{dp} \left\{ \frac {\sinh {A \sqrt p}}{\sinh{B \sqrt p}} \right\} + \left\{ \frac {\sinh {A \sqrt p}}{\sinh{B \sqrt p}} \right\} \right]_{p=0}[/tex]

From WolframAlpha I evaluated:

[tex] \frac{d}{dp} \left\{ \frac {\sinh {A \sqrt p}}{\sinh {B \sqrt p}} \right\} = \frac {1}{2 \sqrt p \sinh{B \sqrt p}} \left\{A \cosh {A \sqrt p} – B \sinh {A \sqrt p} \frac {\cosh {B \sqrt p}}{\sinh {B \sqrt p}} \right\} [/tex]

So the term for s = 1 becomes

[tex]kt \left[ \frac {\sqrt p}{2 \sinh{B \sqrt p}} \left\{A \cosh {A \sqrt p} – B \sinh {A \sqrt p} \frac {\cosh {B \sqrt p}}{\sinh {B \sqrt p}} \right\} + \left\{ \frac {\sinh {A \sqrt p}}{\sinh{B \sqrt p}} \right\} \right]_{p=0}[/tex]

Now I also provide that:

[tex]\sinh{z} = z \prod_{k=1}^\infty \left( 1+\frac{z^2}{k^2 \pi^2} \right)[/tex] and [tex]\cosh{z} = \prod_{k=1}^\infty \left( 1+\frac{4z^2}{(2k-1)^2 \pi^2} \right)[/tex]

So it follows that:

[tex]\sinh{Q \sqrt p} = Q \sqrt p \prod_{k=1}^\infty \left( 1+\frac{Q^2 p}{k^2 \pi^2} \right)[/tex] and [tex]\cosh{Q \sqrt p} = \prod_{k=1}^\infty \left( 1+\frac{4Q^2 p}{(2k-1)^2 \pi^2} \right)[/tex]

Where Q = A or B. Again to save myself the trouble of writing those terms out a million times, I make the substitution:

[tex] \prod_{k=1}^\infty \left( 1+\frac{Q^2 p}{k^2 \pi^2} \right) = \Pi_{Qs}[/tex] and [tex] \prod_{k=1}^\infty \left( 1+\frac{4Q^2 p}{(2k-1)^2 \pi^2} \right) = \Pi_{Qc} [/tex]

Which means:

[tex]\sinh{Q \sqrt p} = \Pi_{Qs} Q \sqrt p [/tex] and [tex]\cosh{Q \sqrt p} = \Pi_{Qc} [/tex]

Alright, making substitutions as described, I am left with the s = 1 term as:

[tex]kt \left[ \frac {1}{2B \Pi_{Bs}} \left\{A \Pi_{Ac} – \frac {A \Pi_{Bc}}{\Pi_{Bs}} \right\} + \frac {A \Pi_{As}}{B \Pi_{Bs}} \right]_{p=0}[/tex]

When p = 0,

[tex] \Pi_{As} = \Pi_{Ac} = \Pi_{Bs} = \Pi_{Bc} = 1 [/tex]

So therefore the S = 1 term reduces to the first term in my solution in my earlier post:

[tex]kt \frac {A}{B} = kt \frac{L-x}{L} [/tex]

Here's where things start to go wrong:

I take a similar approach to the s = 2 term, which is (making the same substitution of A and B as above, and using the rule for derivative of products of functions):

[tex]\frac {k}{2} \left[ \frac{d^2}{dp^2} \left\{ \frac {p \sinh {A \sqrt p}}{\sinh {B \sqrt p}} \right\} \right]_{p=0} = \frac {k}{2} \left[2 \frac{d}{dp} \left\{ \frac { \sinh {A \sqrt p}}{\sinh {B \sqrt p}} \right\}+ p \frac{d^2}{dp^2} \left\{ \frac { \sinh {A \sqrt p}}{\sinh {B \sqrt p}} \right\} \right]_{p=0} [/tex]

The first derivative part is identical to above and the second derivative part I got from WolframAlpha as:

[tex] \frac{d^2}{dp^2} \left\{ \frac { \sinh {A \sqrt p}}{\sinh {B \sqrt p}} \right\}=\frac{1}{\sinh{B \sqrt p}} \left( \frac {A^2 \sinh{A \sqrt p}}{4p}-\frac{A \cosh{A \sqrt p}}{4p^{3/2}} \right)+\sinh{A \sqrt p} \left( \frac{B^2}{4p \sinh^3{B \sqrt p}}+\frac{B^2 \cosh^2{B \sqrt p}}{4p \sinh^3{B \sqrt p}}+ \frac {B \cosh {B \sqrt p}}{4 p^{3/2} \sinh^2{B \sqrt p}} \right) - \frac {AB \cosh{A \sqrt p} \cosh {B \sqrt p}}{2p \sinh^2{B \sqrt p}}[/tex]

Therefore, the entire s = 2 term is:

[tex] \frac {k}{2} \left[ \frac {1}{\sqrt p \sinh{B \sqrt p}} \left(A \cosh {A \sqrt p} – B \sinh {A \sqrt p} \frac {\cosh {B \sqrt p}}{\sinh {B \sqrt p}} \right)+\frac{p}{\sinh{B \sqrt p}} \left( \frac {A^2 \sinh{A \sqrt p}}{4p}-\frac{A \cosh{A \sqrt p}}{4p^{3/2}} \right)+p \sinh{A \sqrt p} \left( \frac{B^2}{4p \sinh^3{B \sqrt p}}+\frac{B^2 \cosh^2{B \sqrt p}}{4p \sinh^3{B \sqrt p}}+ \frac {B \cosh {B \sqrt p}}{4 p^{3/2} \sinh^2{B \sqrt p}} \right) - \frac {AB \cosh{A \sqrt p} \cosh {B \sqrt p}}{2 \sinh^2{B \sqrt p}} \right]_{p=0} [/tex]

Replacing the hyperbolic functions with the infinite product expansions, I get:

[tex] \frac {k}{2} \left[ \frac {1}{Bp \Pi_{Bs}} \left(A \Pi_{Ac} – \frac {A \Pi_{As} \Pi_{Bc}}{\Pi_{Bs}} \right) + \frac{\sqrt p}{B \Pi_{Bs}} \left( \frac {A^3 \Pi_{As}}{4 \sqrt p}-\frac{A \Pi_{Ac} }{4p^{3/2}} \right) + p^{3/2} A \Pi_{As} \left( \frac{B^2}{4p (B \sqrt p)^3 \Pi_{Bs}^3} + \frac{B^2 \Pi_{Bc}^2 }{4p (B \sqrt p)^3 \Pi_{Bs}^3} + \frac {B \Pi_{Bc}}{4 p^{3/2} (B \sqrt p)^2 \Pi_{Bs}^2} \right) - \frac {AB \Pi_{Ac} \Pi_{Bc} }{2 (B \sqrt p)^2 \Pi_{Bs}^2} \right]_{p=0} [/tex]

Which I simplify to:

[tex] \frac {k}{2} \left[ \frac {1}{Bp \Pi_{Bs}} \left(A \Pi_{Ac} – \frac {A \Pi_{As} \Pi_{Bc}}{\Pi_{Bs}} \right) + \frac {A^3 \Pi_{As}}{4B} - \frac{A \Pi_{Ac}}{4Bp \Pi_{Bs}}+\frac{A \Pi_{As}}{4Bp \Pi_{Bs}^3}+\frac{A \Pi_{Bc}^2 \Pi_{As}^2}{4Bp \Pi_{Bs}^3}+\frac{A \Pi_{Bc} \Pi_{As}}{4Bp \Pi_{Bs}^2}-\frac{A \Pi_{Ac} \Pi_{Bc}}{2Bp \Pi_{Bs}^2} \right]_{p=0} [/tex]

Well, this is where I get stuck. When p = 0, all the infinite product terms go to 1, and I would left with:

[tex] \frac {k}{2} \left[ \frac { \left(A – A \right)}{Bp} + \frac {A^3}{4B} - \frac{A}{4Bp }+\frac{A }{4Bp}+\frac{A}{4Bp}+\frac{A}{4Bp}-\frac{A}{2Bp} \right]_{p=0} [/tex]

In which case all the terms with p in the denominator would cancel out and I’m left with:

[tex] \frac {k}{2} \frac {A^3}{4B} = \frac {k (L-x)^3}{8LD}[/tex]

But that’s obviously not the correct value for the s = 2 term because it doesn’t satisfy the boundary conditions. Also, that conveniently ignores the fact that I had a bunch of p terms in denominators before I canceled them out when setting the infinite product factors to 1. Which can't be right.

Well, that’s my working. If anyone can help me figure out where I went wrong, I’d be grateful.

 

[Fred Wright]

I admire your Herculean effort! I could be wrong but I get the subsidiary solution $$
C=\frac{k}{s^2}\left [1-\frac{\sinh(\sqrt(\frac{s}{D})x)}{\sinh(\sqrt(\frac{s}{D})L)}\right ]
$$

 

[Corribus]

Hi Fred. Thanks for taking a look. I don't think you're wrong - but (assuming your s is my p) this seems to be equivalent to my solution, which can be shown by substituting exponential functions for the sinh functions.