The diffusion equation with time-dependent boundary condition

In summary, the conversation discusses the use of Laplace Transform method to solve a 1D diffusion problem with a time-dependent boundary condition. The method involves solving the subsidiary equation and applying boundary conditions to arrive at the solution. However, there seems to be an error in the inverse transform method used, leading to a wrong solution. Suggestions are requested to identify the error. The conversation also touches upon the use of the hyperbolic sum formulae and the residue method in solving the problem.
  • #1
10
4
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.
 
  • Like
Likes Delta2
Physics news on Phys.org
  • #2
Hi. I am curious to know how you get sinh formula by Laplace transformation.
 
  • #3
anuttarasammyak said:
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.
 
  • Like
Likes Delta2 and anuttarasammyak
  • #4
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].

Corribus said:
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]
 
  • #5
Corribus said:
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.
 
Last edited:
  • #6
@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.
 
  • #7
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.
 
  • Like
Likes Delta2
  • #8
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 ]
$$
 
  • #9
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.
 

1. What is the diffusion equation with time-dependent boundary condition?

The diffusion equation with time-dependent boundary condition is a partial differential equation that describes the behavior of a diffusing substance over time, taking into account changes in the boundary conditions (such as concentration or temperature) over time.

2. What is the physical significance of the diffusion equation with time-dependent boundary condition?

The diffusion equation with time-dependent boundary condition is used to model a wide range of physical phenomena, such as the spread of pollutants in the environment, the diffusion of heat in a material, and the movement of chemicals in biological systems.

3. How is the diffusion equation with time-dependent boundary condition solved?

The diffusion equation with time-dependent boundary condition can be solved using various numerical methods, such as finite difference, finite element, or spectral methods. These methods discretize the equation and solve it iteratively to approximate the solution.

4. What are some applications of the diffusion equation with time-dependent boundary condition?

The diffusion equation with time-dependent boundary condition has numerous applications in various fields, including chemistry, physics, biology, and environmental science. It is used to study diffusion processes in gases, liquids, and solids, and to model the behavior of complex systems.

5. What are the limitations of the diffusion equation with time-dependent boundary condition?

The diffusion equation with time-dependent boundary condition assumes that the diffusing substance is well-mixed and that the boundary conditions are constant over time. This may not accurately reflect real-world situations, and more complex models may be needed to capture the full behavior of the system.

Similar threads

Replies
13
Views
2K
Replies
2
Views
1K
Replies
4
Views
998
  • Differential Equations
Replies
3
Views
2K
  • Differential Equations
Replies
22
Views
2K
  • Differential Equations
Replies
1
Views
1K
  • Differential Equations
Replies
1
Views
2K
  • Differential Equations
Replies
8
Views
3K
  • Differential Equations
Replies
2
Views
1K
Replies
1
Views
1K
Back
Top